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Abstract 


Hartung, Lin C. Nonequilibrium Radiative Heating Prediction Method for 
Aeroassist Flowfields with Coupling to Flowfield Solvers (Under the direction 
of Dr. F. R. DeJarnette) 


A method for predicting radiation absorption and emission coefficients in ther- 
mochemical nonequilibrium flows is developed. The method is called LORAN: the 
Langley Optimized RAdiative Nonequilibrium code. It applies the smeared band 
approximation for molecular radiation to produce moderately detailed results and is 
intended to fill the gap between detailed but costly prediction methods, and very 
fast but highly approximate methods. The optimization of the method to provide 
efficient solutions allowing coupling to flowfield solvers is discussed. Representative 
results are obtained and compared to previous nonequilibrium radiation methods, as 
well as to ground- and flight- measured data. Reasonable agreement is found in all 
cases. A multi-dimensional radiative transport method is also developed for axisym- 
metric flows. Its predictions for wall radiative flux are 20 to 25 percent lower than 
those of the tangent slab transport method, as expected, though additional investiga- 
tion of the symmetry and outflow boundary conditions is indicated. The method has 
been applied to the peak heating condition of the Aeroassist Flight Experiment ( AFE) 
trajectory, with results comparable to predictions from other methods. The LORAN 
method has also been applied in conjunction with the computational fluid dynamics 
(CFD) code LAURA to study the sensitivity of the radiative heating prediction to 
various models used in nonequilibrium CFD. This study suggests that radiation mea- 
surements can provide diagnostic information about the detailed processes occurring 
in a nonequilibrium flowfield because radiation phenomena are very sensitive to these 


processes. 
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Summary 


iii 


When a spacecraft enters the atmosphere of a planet, it pushes some of the gas 
particles it encounters with it. As the atmosphere becomes thicker, the gas particles in 
this “shock layer” begin to collide with each other. This converts some of the energy 
of motion that these particles gained from the spacecraft into internal motions of the 
electrons and atoms inside the gas particles. When the gas density is high enough, 
this process reaches an equilibrium state and can be simply described. When the 
gas density is not so high, the transfer of energy proceeds slowly and remains out of 

equilibrium. 

The internal changes in electron and atom energy can be accompanied by a photon 
of radiative energy being emitted or absorbed by a gas particle. Predicting the number 
and energy of these photons in a nonequilibrium gas is the problem to be solved here. 
These photons travel through the gas at the speed of light and may end up depositing 
their energy on the spacecraft surface, thus adding to the heating which exists because 
of gas friction. If enough energy is deposited the surface may begin to melt or burn. 
Heating predictions allow the spacecraft heatshield to be designed accordingly. 
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1 Introduction 

A significant portion of the heating experienced by a blunt spacecraft encounter- 
ing a planetary atmosphere at high speed can be due to radiation. The radiative 
heating can therefore be a strong driver on the design of a spacecraft heatshield. In 
the past, most spacecraft encountering planetary atmospheres were on entry trajecto- 
ries descending rapidly into the dense lower layers of the atmosphere. The heat load 
they experienced could therefore be accurately calculated using equilibrium methods. 
Mission scenarios currently under study, such as Aeroassist Space Transfer Vehicles 
(ASTV) designed to rendez-vous with the International Space Station Freedom, and 
aerobrakes used for orbit insertion in planetary exploration, present a different situa- 
tion. Walberg [1] provides a summary of such missions. These vehicles use the upper 
atmosphere to obtain the necessary velocity change for orbit transfer and may spend 
significant amounts of time in low density regions where nonequilibrium conditions 
prevail. 

Numerous Computational Fluid Dynamics (CFD) codes are being developed to 
solve the Navier Stokes equations for this nonequilibrium regime. These codes predict 
the convective heating rates for AST Vs. A few studies of the nonequilibrium radiative 
heating problem have also been made [2, 3, 4]. Nelson [5] discusses some ot the 
problems and preliminary results from nonequilibrium radiative heating analyses. 
However there remains a need for a relatively fast and accurate method which can 



be coupled to a flowfield solution in cases of heavy radiation, or used uncoupled m 
cases with less significant radiation. There has also been some evidence in recent 
years [4] that the one-dimensional tangent slab approximation, which is standard in 
computing radiative transport, may not be adequate for some of the vehicle shapes 
and applications under consideration. It will certainly not suffice to compute radiative 
transport in the wake region of such vehicles where payloads will generally be located. 

The objective of the present work is to develop an accurate, efficient method for 
calculating radiative heat transfer under nonequilibrium conditions, to initiate the 
development of a multi dimensional transport algorithm for such radiative transfer, 
and to present some initial results. 

To develop a method for calculating nonequilibrium radiation, the necessary radi- 
ation property models were first derived for a nonequilibrium gas. In nonequilibrium, 
expressions for both absorption and emission are required since the equilibrium rela- 
tionship between the two does not apply. This theoretical development is presented in 
Ch. 3. These nonequilibrium radiation models were then implemented in a comp 
code dubbed LORAN: the Langley Optimized RAdiative Nonequilibrium code. As its 
name suggests, this code is optimized to provide radiation predictions efficiently, by 
a judicious selection of spectral points determined from the distribution of radiative 
transitions. Parametric studies were then performed to determine the minimum set 
of optimized spectral points which would provide accurate radiative heat flux predic- 
tions. This process and its results are described in Ch. 5. The optimized method 
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is compared to flight- and ground-based measurements, and to other nonequilibrium 
radiation methods, in Ch. 6. The agreement is generally good. 

To compute radiation transport without applying the tangent slab (1-D) approx- 
imation, a full 3-D modified differential approximation (MDA) was adopted. The 
derivation of this method in nonequilibrium is summarized in Ch. 4, while the nu- 
merical details of its implementation are presented in App. B. This method is cur- 
rently limited to axisymmetric flowfields, though the extension to 3-D flows should 
be straightforward. The method has been demonstrated for two nonequilibrium flight 
conditions and compared to results from the tangent slab approximation in Ch. 6. 
The predicted wall radiative heating is 20 to 25 percent lower than the tangent slab 
result, in line with what would be expected. There are, however, indications that the 
boundary conditions of the MDA method can be improved. The radiati\e coupling 
term, V9R» varies much more smoothly in the MDA solution than in the tangent slab 
prediction, suggesting that the MDA method may enhance the stability of coupled 

solutions. 

Chapter 7 applies the LORAN radiation model to a study of the sensitivity of 
radiative heating predictions. Semi-empirical models for the rate controlling tempei- 
ature in dissociation reactions, the cross section for vibrational to translational eneig\ 
exchange, and the amount of energy exchange in dissociation, as well as various chem- 
ical kinetics models used in nonequilibrium CFD were varied and the effect on the 
predicted radiative heating levels was examined. The results show' that radiation is 
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very sensitive to these models through their influence on the vibrational temperature. 
Studies of radiative heating may therefore provide additional insight into the valid- 
ity of the semi-empirical models which have been developed for nonequilibrium CFD 

codes. 

The nonequilibrium radiation and transport methods developed for this work pre- 
dict radiative heating to aeroassist vehicles with accuracy similar to that of detarled 
codes. The computer time and storage requirements are significantly reduced, how- 
ever. The LORAN method is therefore suitable for use in parametric studies and for 
coupled solutions with nonequilibrium Naviet Stokes codes. Further work on these 
models is expected to speed up solutions still more. 
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2 Historical Review 

2.1 Topical and Literature Survey 

The previous work which has relevance to this problem can be subdivided into 
several subject areas: nonequilibrium absorption and emission in individual radiating 
mechanisms, nonequilibrium excitation processes, and transport models. 

2.1.1 Nonequilibrium Absorption and Emission 

Emission and absorption of radiation in a gas result from transitions between 
energy levels. In atomic species, all electronic energy levels as well as free electrons 
may be involved in radiative transitions. For molecules, radiation also occurs due to 
transitions between rotational and vibrational energy levels. Atomic radiation arises 
from bound-bound transitions (atomic lines), bound-free transitions (photoionization 
or radiative deionization), and free-free transitions (Bremsstrahlung radiation). Here 
“bound” and “free” refer to the state of the electrons involved in the radiation- 
inducing transition. Molecular radiation is more complex because of the large number 
of energy levels available for transitions. The resulting structure of closely-spaced lines 

is referred to as a molecular band system. 

The complete spectral variation of emission and absorption is contained in the 
absorption coefficient, k„, and the emission coefficient,;'. The absorption coefficient 
is obtained from the product of the population, iV„, of an energy level and its radiative 


6 


absorption cross section. <T„„. summed over all the energy levels available: 

energy levels n 

The emission coefficient is found from a similar expression, in which the amount of 
energy emitted, fuw also appears, and the absorption cross sect, on is replaced by 

the transition probability for emission, A nn '- 

J e = E Nn Z hu nn -A nn > C 2 - 2 ) 

energy levels n n ' <n 

In equ, librium, the electronic energy level populations are determined as a func- 
tion of temperature according to a Boltzmann distribution. The rotational and vi- 
brational energy levels of each electronic state of a molecule are populated according 
to the equilibrium partition function. Species concentrations in equilibrium are re- 
lated by Saha’s equation. These equilibrium relations allow radiative properties to 
be referenced to temperature and the concentrations of a few major species, allowing 
the development of simple radiation step models. Also, under equilibrium conditions 
the absorption and emission coefficient are related according to Kirchhoff s law. so 
that Eq. 2.2 is not needed. 

In nonequilibrium, these simple relations for the energy state populates and 
species concentrations no longer apply. The problem is now to determine the nonequi- 
librium energy level populations. The rotational and vibrational states in a molecule 
may be assumed to be populated according to distinct rotational and vibrational tem- 
peratures, respectively, if the nonequilibrium is not too severe. The electronic states. 
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however, equilibrate more slowly so that a single electronic temperature cannot be 
defined. In this situation, the energy level populations in the absorption and emis- 
sion coefficients must be calculated individually at every point where they are needed 
according to the specific local conditions. In addition, the emission coefficient can- 
not be calculated by Kirchhoff’s law. The radiative cross sections, however, depend 
only on the configuration of an individual atom or molecule. They are the same for 
equilibrium and nonequilibrium conditions. 

A number of sources present expressions for the equilibrium absorption and/oi 
emission coefficients of the individual radiating mechanisms. For example, Zel dovich 
and Raizer [6], provide a comprehensive classical development with discussion and 
references to quantum corrections. Other classic texts present similar results [7]. Ai- 
maly [8] provides practical expressions for the continuum absorption coefficients of 
atoms and ions, following the method of Biberman and Norman [9]; which itself fol- 
lows their earlier work [10] and that with Ulyanov [11]. The latter work is an attempt 
to correct the classical hydrogenic model of the atom to account for the actual con- 
figuration of a radiating species and its influence on the absorption coefficient. Jones 
et al. [12] provide a critical review of many of the important radiating mechanisms 

and species. 

A number of other researchers have also considered the contributions of various 
parts of the spectrum in a more qualitative sense for conditions of interest in re-entr\ 
studies. Horton [13] has studied the importance of the ultraviolet portion of the 



spectrum as a contributor to the radiative heating of re-entry bodies. He has also 
studied the potential impact of nonequilibrium on the radiative heating [14]. 

Only two approaches to absorption and emission in nonequilibrium have been 
identified. The correction factor method, developed by Clarke and Ferrari [15] and 
Chapin [2], and now in use by L. A. Carlson and others [16, 17, 18], incorporates a 
simple equilibrium step model for absorption with correction factors to account for 
the nonequilibrium of the excited state populations. The detailed line by line method 
developed by Whiting et al. [19] and implemented for nonequilibrium only in Park s 
NEQAIR code [20] treats individual rotational lines in the molecular band spectrum 
as well as individual atomic lines and atomic continuum processes. In this method the 
nonequilibrium absorption coefficient is obtained by solving for the nonequilibrium 
populations of the energy levels using the quasi-steady-state (QSS) method described 
in the next section. The nonequilibrium emission is obtained by calculating a quasi- 
equilibrium or excitation temperature for each radiative process, allowing the use of 
Kirchhoff’s law. 

2.1.2 Nonequilibrium Excitation Processes 

To calculate the nonequilibrium populations appearing in the absorption coeffi- 
cient, the rates of population and depopulation of energy levels must be known or 
assumptions must be made about how the states are populated. Excitation rate data 
are not widely available because of the experimental difficulties involved in obtain- 
ing them. In the future, computational chemistry may be applied to predict these 
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rates, but such calculations are costly and therefore computer- and time-limited. Park 
[3, 20] has collected a set of rates for the principal transitions involved in air radiation 
and incorporated them into a quasi-steady-state (QSS) model of electronic excitation. 
The QSS model assumes that the rates of population and depopulation of an energy 
level are much larger than the net rate of change of the level's population. This allows 
the excitation calculation to be uncoupled from the flowfield solution. In flow regions 
with strong gradients (such as inside captured shocks and in boundary layers) this 
approximation loses accuracy because it neglects the excitation history. 


2.1.3 Radiative Transport 

The problem of radiative transport is well-known. It requires solving an integro- 
differential equation for each frequency and for each direction in the flowfield. This 
complexity arises because unlike a temperature field in conduction problems, the ra- 
diation intensity in any direction is independent of that in any other direction: and 
that at a given frequency (energy) is separate from that at any other frequency. The 
complete equation of radiative transport has been derived in many sources, including 
ZePdovich and Raizer [6], Ozisik [21], and Vincenti and Kruger [22]. Using Ozisik s 
notation, the governing equation for the spectral radiative intensity, /*,, in a paitici 

pating medium is 


i dhjs, n, t) 

c dt 


a/„(3, n,o 

ds 






(2.3) 
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where k' u is the absorption coefficient corrected for induced emission. 

The first term on the left-hand-side is immediately omitted in most radiative heat 
transfer studies because the speed of light, c, is large and the time derivative of 
intensity is small in most applications of interest. The scattering terms (containing 
7 * and p) can also be neglected, since the emphasis of the present study is on non- 
ablating, low-density flows. In such flows there is almost nothing to induce radiation 
scattering, and 7 * « 0. With these two simplifications, the transport equation reduces 
to: 

dL(s, n. 1 ) + n, t) = ms, t) (2.4) 

OS 

where t is now a parameter. This equation can be solved numerically with the as- 
sumption of a one-dimensional medium. For multi-dimensional radiation transport, 
additional simplifications will be necessary. The most common radiative transport 
models used to date are the optically thin model and the tangent slab model. The 
first assumes a transparent gas, so that the transport problem reduces to a simple 
summation across the region of interest. The second assumes that all properties vary 
in only one direction, so that the radiative transport equation is reduced to a one- 
dimensional integro-differential equation. A typical discussion of the tangent slab 
method can be found in Ozisik [21] or in Nicolet [23, 24]. An extension of Nicolet's 
method is given by Bolz [25]. 

Several other approximate methods for the solution of the radiative transport 
problem have been developed. The optically thick approximation assumes that emit- 
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ted radiation is reabsorbed within a very small distance from the point of emission 
[21. Sec. 9.2]. This is not the case for blunt body flowfields, except perhaps at the 
centers of very strong atomic line transitions, so this method is not applicable to the 
current problem. Another series of methods, known variously as moment methods, 
discrete ordinate methods, or spherical harmonics methods, obtain solutions by ap- 
proximating the radiative intensity with a series of functions of increasing order [21. 

Sec. 9.8]. 

Multi-dimensional solutions to the transport problem have been developed by 
Truelove [26] and Modest [27], but these assume equilibrium flow and the validity 
of Kirchhoff’s law. In addition, Truelove assumes a gray gas, in which the absorp- 
tion does not vary with wavelength. Cheng and Ozisik [28] have developed a three- 
dimensional algorithm for stagnation flows, but it uses the diffusion approximation tor 
radiative transfer which is valid only in highly absorbing flows. No multi-dimensional 
transport algorithm currently in use for flowfield solutions has been identified, though 
an effort toward this end is underway by Edwards et al. [29]. 

Equally important to the development of radiative transport algorithms may be 
the formulation of rules governing the intervals at which radiation must be computed. 
Especially in a coupled solution, it is preferable to minimize the number of times the 
costly radiation model is invoked. Bolz [30] has developed an algorithm for selecting 
a radiation transport subgrid in an equilibrium flowfield by evaluating a weighted av- 
erage of the derivatives of important flow variables. This avoids calculating radiation 
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at every grid point in a flowfield region where all the properties are changing very 
slowly. 

2.2 Computer Code Survey 

A large number of computer codes for radiation prediction have been developed 
over the years. Unfortunately, many of these codes fell into disuse during the L970's 
when hypersonic research was limited. As a result, most are no longer available for 
current study without extensive redevelopment. This survey, therefore, will consider 
only those codes known to be available with minimal reconstruction required. 

2.2.1 Equilibrium Radiation Codes 

Gray gas 

Gray gas models are the simplest of all radiation models, assuming that emission 
and absorption do not vary with frequency. This assumption is violated for blunt 
body flowfields, so these models have been rejected as candidates for the current 

effort and will not be discussed further. 

Step models 

A number of researchers have developed step models, in which the absorption 
coefficient is considered constant in several wavelength intervals. Olstad [31] has 
developed an eight step model, while Zoby et al. [32] used a 58-step model to study 
outer planet entry. These step models depend on the assumption of equilibrium. 
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They make extensive use of curve-fits and refer radiation mechanisms to a ground 

state in order to obtain the desired simplicity. 

The eight-step model developed by Olstad [31] is still available. This model makes 

use of curve-fits for the radiative cross sections assuming equilibrium conditions. It 
requires input of the concentrations of five species (N, O, N„ 0, and and uses 
these to infer the concentrations of minor species, again assuming eqmlibtium. The 

details of this model can be found in the reference. 

Step models for radiation emphasise simplicity and speed at the expense of ac 
curacy. This trade-off is often acceptable for engineering methods which provide 
estimates of the radiative heat load. These simple models are ill-suited to nonequi- 
librium, however, where properties are no longer simple functions of temperature 
and pressure but depend also on the concentrations of the various species. In .ha, 
situation, the simple relationships used to develop step models no longer hold. 


RAD/EQUIL 

The RAD/EQUIL code developed by Nicolet [23, 33) is still in use in its original 
form. It also serves as the basis for at least two coupled equilibrium invisc.d flow 
codes: RIFSP, which computes a stagnation line for a blunt body, and RAIF, uhi 
computes a complete blunt forebody [34]. This method is a hybrid between the step 
models and detailed models. It computes continuum processes from curve-fits for 
equilibrium, and reports the results a. a finite number of points (generally 20-30). 
The line radiation is computed in greater detail, with each line considered being 
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represented by fifteen spectral points. However, these detailed results are integrated 
by the code, and presented at a small number of intervals (detailed output can also 
be obtained, but is not the default). This process provides results of good accuracy 
and allows concise reporting. 

2.2.2 Nonequilibrium Radiation Codes 
RADMC 

This nonequilibrium correction factor method is currently being developed by a 
number of researchers [16, 17, 18, 35]. It will be referred to collectively as the RADMC 
method. In this method, the eight-step model for absorption developed by Olstad 
[31] is adopted to describe the radiation spectrum. Approximate correction factors 
for nonequilibrium energy level populations are then developed. These depend on 
modeling the excitation with a small number of energy levels. The result is a very 
fast but approximate nonequilibrium radiation model. 

NEQAIR 

The detailed line-by-line approach developed by Whiting et al. [19] is implemented 
for nonequilibrium in the NEQAIR code of Park [20]. This method includes detailed 
descriptions of the radiative transitions, computing each individually. The nonequilib- 
rium electronic state populations are calculated using the QSS approximation. This 
is a very detailed but computationally costly model for nonequilibrium radiation. 

on the order of 10 5 to 10 6 spectral points [29, 36]. 


Typical calculations require 
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Direct Simulation Monte Carlo 

An emerging method models radiation in a Direct Simulation Monte Carlo (DSMC) 
flowfield solution. The radiation modeling for DSMC is still being refined but has 
great potential, particularly for nonequilibrium low density flows. The power of 
DSMC lies in its direct modeling of the physical processes involved, and therefore 
suffers most from uncertainties in what these models and their parameters should be. 

It also requires large amounts of computer time, though advancements are reducing 
this considerably. The computer time required is proportional to the flow densitt. 
however, restricting the usefulness of DSMC for lower altitudes. Moss and Price give 
some initial results for the Aeroassist Flight Experiment (AFE) vehicle in both fore- 
body and wake regions [37] using a method developed by Bird [38], A. B. Carlson 
[39, 40] has also studied a shock tube condition with a DSMC method including 

radiation. 

2.3 Survey of Experimental Results 

2.3.1 Ground-Based Data 

The principal source of ground-based data is a series of shock tube experiments 
performed by the Avco-Everett Research Laboratory around the early 1960's. These 
results are reported in a number of sources, including Teare et al„ [41], Alien et al. 
[42, 43. 44, 45], and Kivel [46], The conditions tested include velocities between 6 
and 10 km/sec. Only the high end of this range is of interest for the current problem. 
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Experiments were also performed by Page [47, 48, 49] and Nerem et al. [50. ol. 
52, 53, 54, 55, 56, 57, 58, 59, 60]. These shock tube results are difficult to use because 
the data were reinterpreted in the years following the experiments. No attempt will 
therefore be made to match this data here. 

2.3.2 Flight Data 

Project FIRE 

The FIRE flight project of the mid- 1960’s was conducted to determine the heat 
load on an Apollo-type vehicle entering Earth’s atmosphere at 11.4 km/sec [61, 62]. 
It is one of the only sources of flight radiative heating data applicable to the present 
area of study and will therefore be emphasized in verifying the development of this 
radiative heating prediction method. 

Future Flight Projects 

Looking to the future, NASA’s Aeroassist Flight Experiment is anticipated to 
provide radiative heating data for a blunt body on an aerobraking trajectory in the 
earth’s atmosphere. The heat shield design is non-ablating, and both the forebody 
and wake regions are planned to be instrumented. The data from this project should 
provide a valuable verification set for radiation prediction methods. 


3 Radiation Theory 


Absorption and emission coefficients describing the radiation phenomena occur- 
ring in a blunt body flowfield can be developed under conditions of thermochemical 
nonequilibrium for each participating transition between energy levels. The develop- 
ment for the absorption coefficient corresponding to each type of transition is pre- 
sented below. A similar development for the emission coefficients is given in Sec. 3.2. 

3.1 Absorption 

3.1.1 Atomic Mechanisms 

Bound-Bound Transitions 

The model for radiation absorption in bound-bound transitions used in the present 
study is that developed by Nicolet [23] for the RAD/EQUIL code. This model resolves 
a line by distributing about 15 spectral points starting from the line center as shown 
in Fig. 3.1. It has been modified as required for nonequilibrium. Following Nicolet ‘s 
development, the absorption coefficient in a line is given by 

*„ = (3.1) 

me 

where f tu is the oscillator strength of the transition, Ni is the nonequilibrium popu- 
lation of the lower state of the transition, and b is the line shape function which may 
depend on a large number of variables. The constants e, m, and c are the electron 
charge, electron mass, and speed of light, respectively. 
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Figure 3.1. Distribution of Spectral Points for an Atomic Line 

The line shape is a normalized function which can be described by one of a numbei 
of theoretical models depending on the predominant broadening mechanism. Two 
such line shapes are considered here. The Lorentz shape describes lines in which Staik 
broadening by electron impact dominates, as is the case for heavy atomic species. This 
line shape is given as 

u _ ( 3 . 2 ) 

{v - vcl) 2 + ( 7 s ) 2 

where the dependence of b on the gas conditions has been omitted. Here -7 s is the 
Stark (half) half width, and v C l is the frequency of the line center. The (half) half 
width is the half-width of the line at half its height, and is a common descriptor for 
line shapes. 

The (half) half widths have been calculated by a number of authors, most notably 
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Griem [63]. Page et al. [64] have noted that these line widths vary with a power n of 
temperature T as 

7 s (D« (£)V( r o) ,3 :5) 

thus allowing use of a single input reference at the reference temperature T 0 for 
Stark broadening. In nonequilibrium, the temperature in this equation should be the 
electron temperature T e . 

The mechanism of resonance broadening (also called collision or pressure broaden- 
ing) yields Lorentz line shapes as well. This mechanism may exceed Stark broadening 
for situations with low ionization. Nicolet [23] gives the approximate (half) half width 


as 


/v\ 

e 2fR 1 

UJ 

2innv R 


N a 


(3.4) 


where N a is the number of perturbing atoms per unit volume, and gi and g u are the 
degeneracies of the upper and lower states, respectively. The effective Lorentz line 

width is then the sum of 7 s and 7 *. 

The second type of line shape is the Doppler profile, which describes line broad- 
ening due to thermal motion of the atoms. This Gaussian line shape is given as 


by — D 


1 In 2 

(u — ucl) 2 In 2 

/ ex P 

j 7T 

r (7 d ) 2 J 


[3.51 


The Doppler (half) half width is given by 


D _ ''Cl 1 -kTt In 2 
C I m s 


(3.6) 
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where k is Boltzmann's constant and m, is the mass of the radiating part.de of spec.es 
s. The temperature is the heavy particle translational temperature, since the effect 

is due to the thermal motion of the atoms. 

The Lorentz and Doppler line shapes can be combined into a single curve known 

as the Voigt line shape which requires significantly more computation. Under the 
conditions of interest in this study, the two (half) half-widths have been found to be 
of different orders of magnitude. Therefore, it is felt that the use of one or the other 
of the simpler line shapes should be adequate. Once the choice is made for each given 
condition, the absorption coefficient in a line can be computed from Eq. 3.1, using 
input data for the oscillator strength and the line width and predicted nonequilibrium 

energy level populations. 

Bound-Free Transitions 

A model for bound-free transitions can be developed assuming a hydrogen-like 

atom. A good source of such a development for an equilibrium gas is Zel'dovich 

and Raizer (6, Sec.5.5], The radiative cross sections developed there apply also to 

nonequilibrium since the configuration of an individual atom is unchanged. From 

that reference therefore the cross section for an atomic energy level n is. 

1 647 r 4 mZ 4 e 10 ( 3 , 7 ) 

<7 ‘ /n _ \/3 3 v 3 ch 6 n 5 

where Z is the charge number (Z = 1 for neutral species) and h is Planck’s constant. 
Once the cross section is known it is a simple matter to compute the spectral ab- 


sorption coefficient for the bound-free process, given the nonequilibrium populations 
of the electronic energy levels of the atom. The relation is 

oo 

K„=£]An0Vn (3-S) 

n* 

The lower limit n* indicates the photoionization threshold. This is the lowest atomic 
energy level from which a photon of frequency v can detach an electron. The upper 
limit in practice has a finite value corresponding to ionization of the atom, 

While quite accurate, Eq. 3.7 can still be improved. In particular, a quantum 
correction [6, p. 266] may be applied to it. However, as stated in the reference, 
the impact of this correction is negligible for most cases of practical interest. More 
significant is the error in this expression for low quantum states near the ground 
atomic state. For these states, the hydrogen-like assumption may induce considerable 
error. For many atoms, the cross section of the ground state is well known from 
experimental data. This information can be used for transitions involving the ground 
state instead of the value predicted by Eq. 3.7. 

Free-Free Transitions 

A model for free-free transitions can be developed using similar semi-classical 
methods. Only transitions due to the proximity of an atomic ion are considered. 
Free-free transitions caused by the presence of neutral atoms and molecules are as- 
sumed to be negligible. Zel’dovich and Raizer [6, Ch. 5.3] is again a good source. 
Their result can be rederived without using equilibrium assumptions except for a 
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Maxwell distribution of the energies of the free electrons at the electron translational 
temperature T e . In this case, the only change in the result is that T e replaces T t . The 
equation for the free-free absorption coefficient can therefore be written as 


= 


= i( 


2ir \' 12 Z 


2^6 


3 \3mkT e ) hcmv 3 


N + N e 


( 3 . 9 ) 


where N + and N e are the concentrations of the appropriate ion and of electrons, 
respectively. 

This result could also be modified by a quantum correction, but again the effect 
is minimal for the conditions of interest here. 


3.1.2 Molecular Mechanisms 

The treatment of molecular radiation requires consideration of vibrational and 
rotational transitions within the molecule, as well as electronic transitions; and com- 
binations of the different types of transition. Each such transition contributes a 
discrete line to the radiation spectrum, whose line center frequency is determined by 
the energy differential of the transition. While a line-by-line calculation of molec- 
ular radiation is possible, it does not fit the requirements of the present study for 
a relatively rapid calculation method which can be coupled to a flowfield code. As 
an alternative, therefore, the -‘smeared band” model of molecular radiation will be 
adopted here. In this model, the rotational line structure of the molecular radiation 
is smoothed to provide a continuous variation of the absorption coefficient. Because 
the rotational energy levels are so closely spaced, the lines resulting from rotational 
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transitions are nearly overlapping except in very rarefied conditions. The error intro- 
duced by this approximation should therefore be acceptable under most conditions ot 


interest. 

The basic development of the “smeared band” model can be found in a number oi 
sources. The development given below is an enhancement of that found in Zel do\ ich 
and Raizer [6, Sec. 5.3]. It incorporates higher order expressions for the molecular 
energy levels from Herzberg [65], while eliminating equilibrium assumptions. In what 
follows ZR will refer to the text by Zel’dovich and Raizer. Herzberg will refer to the 

just-quoted reference. 

The energy differentials for the molecular transitions will be established first. For 
electronic transitions, the energy levels are given by £ e , the electron energy of a 
particular energy level measured from the ground state. For vibrational transitions, 
the energy levels are obtained from the fourth order expression given by Herzberg 


E u — he 





+ ^ e y e 






( 3 . 10 ) 


which includes a correction for the anharmonicity of the vibrational energy potential. 
The vibrational levels range from the zero-point energy at u = 0 to some maximum 
r, determined by the dissociation of the molecule. The coefficients in the above ex- 
pression are spectroscopic constants which are available from a number of sources. 
The energy levels for the rotational transitions, including a correction for non-rigid 


oscillators, are 


E r = he \B V J{ J + 1) — D V J (J + 1) ] 


(3.11) 
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where the spectroscopic constant B u has been corrected for rotational effects as 

B v = B e -a e (u + ( 3 - 12) 

and the second rotational spectroscopic constant D v is 


1 4 B e 3 

D v = D e + i3 e (v + -) « D e = — 


(3.13) 

These levels range from zero energy at J = 0 to some maximum J again determined 
by dissociation of the molecule for each vibrational level. 

The complete molecular energy level is 


Etot = E e + E v 4- E r 


(3.14) 


where the individual energy contributions are given by Eqs. 3.10 and 3.11 above. 
Designating the upper level by the superscript « and the lower level by I the wave 
number (inverse wavelength) of a particular transition is 


1 _ El t - E[ ot 
A he 


(3.1o) 


For convenience in manipulating these expressions one special wave number is de- 
fined. This is the wave number with zero rotational quantum number and arbitrary 

vibrational quantum numbers: 


1 = El~Ej ! Et-E[ 

A„u„i he he 


(3.16) 


The formidable expression for the wave number of an individual transition is 
simplified by the selection rules for rotation. These rules are discussed in various 
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sources including ZR and Herzberg, and result in the requirement that A J = J u -J l = 
0. +1, or - 1. In addition, the 0-0 rotational transition is forbidden and in the case 
of transitions between certain molecular states denoted by E all AJ = 0 transitions 
are forbidden. A detailed discussion of these selection rules and of the designation 
of molecular states is outside the scope of the present work. Nevertheless, these 
rules allow the development of three specific wave numbers for the A J = -1. 0, 1 
transitions, which are denoted P, Q, and R respectively. 

-1 = j- + J 1 ' {d[ - d;) + J‘‘ (2 D[ + 2D') 

A/ 3 

+J' 2 (D l e - B[ - £>“ + fl“) + J 1 (-5“ - B[) (3.17) 

— = — b J 1 * {D[ - D*) + f 3 (2 D[ - 2 D*) 

Ag A v u v i 

+J 1 ' (£>; - Bl - d; + b;) + J‘ (b; - el) (3. is ) 

-L = A- + j 1 ' Id 1 , - d;) + J 1 ’ (2 d‘, - a of) + 

A r A v u v i 

f ( D ‘ e _ B[ - 13D“ + B \ “) + J 1 (3B V “ - 12 D\ - B l v ) + 2 - 4 £>“ (3.19) 


Once the molecular band spectra are thus determined, it is necessary to evaluate 


the probability of a transition between levels in order to obtain an expression toi 
the absorption coefficient. This probability is obtained from a quantum mechanical 
analysis. The resulting expression for the transition probability is given for instance 

in ZR: 


A* = A 


Bv u J u 

Av l J l 


64?r 4 3 
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BIiba Qv u v l Pj u j‘ 


(3.20) 
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where q v u L ,i is the Franck-Condon factor for the probability of a vibrational transi- 
tion, and Pjuji is the probability of a rotational transition. D e ^g A is the square of the 
matrix element^ arising from the wave function in the quantum mechanical analysis. 
A is the Einstein coefficient for spontaneous emission, while super- and subscripts A 
and B denote the lower and upper electronic levels, respectively. Using the princi- 
ple of detailed balancing leads to the well-known relationship between the Einstein 
coefficients for emission and absorption (see for instance [66] or ZR) 


Bin = 


Sirhi/^ g t 


(3.21) 


Note that if this relationship is expressed for energy states rather than energy levels 
the degeneracies g u and gi are identically one, by definition. As a practical matter 
the use of energy levels, which may be degenerate, is preferred because it reduces the 
bookkeeping required. 

The absorption coefficient, k/ u , is related to the Einstein coefficient for absorption, 
B iu , by [22] 

k/u = ri[Bi u hi/i u (3.22) 


Substituting Eqs. 3.20 and 3.21 in this last expression yields the desired expression 
for the spectral absorption coefficient: 


gsv^J* 


87T 3 


U Av l J l , Bv u J u 


g Av ijt 3 he 


N A vU‘ B\ iba q v u v i Pjhji F{u) (3.23) 


where F(i/) is a normalized line shape. This expression must be evaluated for each A. 

1 sre2, used by Park in NEQAIR, is the nondimensional value, which is divided by (a 0 e) 2 , the 
squared product of the Bohr radius and the charge of an electron. 



B, u“, u', J u , J‘ in order to obtain the complete absorption spectrum of a molecular 
species. The "smeared band” model averages out the rotational line structure. This 
is done by summing Eq. 3.23 over J u , and introducing an average frequency. 


for the rotational lines. The result is: 


U Av l jl,Bv u 


8tt 3 9b 

3 He q 
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[3.24) 


where use has been made of the relation Ej« pj»jt = 1 which expresses the fact that 
all transitions occur between two levels of the molecule. The degeneracy has been 
expanded into its components: 

9AvJ = 9 a9v9j ( 3 - 25 ) 


where 

9 . = 1 ( 3 . 26 ) 

since there are no degenerate vibrational levels, and 

gj =27+1 (3.27) 

For the vast majority of rotational lines, J 1, and the ratio of rotational degen- 
eracies in Eq. 3.23 is approximately equal to one. This approximation has been used 
even for transitions with small J. While the error introduced in a single transition 
with a low rotational quantum number can be as much as 300 percent, onlv a few 
transitions have low J so the effect on the overall absorption coefficient is small. The 
error incurred by introducing the average frequency Vbv u ,Av 1 f° r rotational lines in 
the v i _ v u band is minimal, on the order of one to seven percent. 
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Proceeding towards the complete absorption coefficient, Eq. 3.24 is now summed 
over the final vibrational level t> u . This yields 


9 B \ r 

TL 

ihcg A 


Here the frequency remains inside the summation and is not averaged, as the variation 
involved can be large. 

The mean value theorem is then used to define an average absorption coefficient 
for a line, by integrating over one line width: 


, on yB A r n2 - - 

v ljl * ' Av l J‘ UelBA 2-1 9v u v‘ l/ Bv u ,Av l ~ k ^ Av i 


where hv is the line width. 

Equation 3.29 is then solved for k„ and summed over the initial rotational and 
vibrational levels. 


3 „ oo oo oo \r 

9b n 2 — Y^ ™ Av l J l 

~ ^elBA 2-* Qv u v l l/ Bv u ,Av t 2-*t \ 

She gA i v u jt 


where the sum term is zero for any transition which does not absorb radiation in the 
frequency interval under consideration. 

Equation 3.30 is the complete formula for the “smeared band” absorption coeffi- 
cient for a molecular species. In this form it contains two approximations: the use of 
an average frequency for the rotational lines in a vibrational band, and the averaging 
over a single line. In order to use this expression in a radiation model, its various 
parts must now be amplified in terms of known quantities. 
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Under conditions of nonequilibrium, the number density in a particular energy 
level is given by the following expression, where the existence of rotational and vibra- 
tional temperatures has been assumed. 


exp (~E v i/kT v ) gjt exp{- Ej,/kT T ) 

• 4 Q V {T V ) Q r {Tr) 


(3.31) 


Q v and Q T are the vibrational and rotational partition functions, respectively. The 
population of the electronic energy level N .4 must be obtained from some nonequi- 
librium excitation calculation. In the present work, Park’s quasi-steadv-state (QSS) 
method [67] has been used. The partition functions are obtained from their definition: 


Q X {T X ) = ^exp {-E x jkT x ) (3.32) 

r=0 

In practice the upper limit of this summation is imposed either by predissociation m 
the rotational case or by numerical limits arising from inaccuracies in the higher oidei 
coefficients of Eq. 3.10 for the vibrational case. The rotational case will be discussed 
in a separate section below. For the vibrational case, the sum may be terminated 
when the contribution of the next term is less than some threshold percentage of the 
partition function. It is also necessary to monitor the size of each term to detect the 
onset of significant inaccuracies in Eq. 3.10. These occur at a high enough vibrational 
quantum number that the sum for the partition function can be truncated at that 
point. 

The next term to be examined in Eq. 3.30 is the At/ term. This term introduces 
another approximation of the “smeared band"’ model. The line width is assumed 
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to be approximately equal to an average line spacing. To evaluate the line spacing, 
consider the three expressions for the line centers, Eqs. 3.17 to 3.19. If the D e terms 
are neglected compared to the B v terms, since they are much smaller, the equations 


reduce to: 


A. 1 + j p (b: - b‘ v ) + j 1 (-b: - b[) 

A p A^U yl 

(3.33) 

i-_ 1 + J‘ 2 (b; - B[) + J‘ (b: - Bl) 

A Q AyUyl 

(3.34) 

i-_ 1 + f (b; - b[) + j 1 (3b; - b[) + 2b; 

A XyUyl 

(3.35) 

For a line with J' » 1, which is the case for most lines, these three expressions 

collapse into the single form: 


i_-~ 1 +j ,2 (b:-b[) 

A c Aw 

(3.36) 

Solving for v and taking the derivative with respect to J, remembering that A J 1 

for adjacent lines, gives 


Ai/ = 2cJ l \B v * — B v i | 

(3.37) 


where the absolute value is taken to ensure that the average spacing is positive. 
Equation 3.36 also yields an approximate expression for J‘ in terms of the wavelength 

or wave number: 

J-l 1-=- (3.3S) 

J - Aw ) B v » - B v i 

Note that since J' 2 > 0, the right-hand-side of this equation must also be positive. 
This implies for B v » > B v i that Aw > A, and for B v - < B v t that Aw < A. These 
requirements are the mathematical statement of the fact that each vibrational band 
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absorbs in only a portion of the spectrum. In evaluating the absorption coefficient for 
a particular band, the wavelength or frequency under consideration must be checked 
against the above conditions to avoid adding superfluous contributions. 

Returning now to Eq. 3.11 for the rotational energy levels, the assumption of 
J' > 1 is again made. This allows the replacements J‘ + 1 J l and 2 J l + 1 ~ 2J ■ 
Equation 3.38 can then be substituted in Eq. 3.11, with the result: 



This expression can then be substituted in the rotational exponent in Eq. 3.31. It 
constitutes the final approximation of the “smeared band model. 

The last term to evaluate in Eq. 3.30 is the averaged frequency, u B v\Av‘- This 
must be evaluated separately for E and non-E transitions according to the rotational 
selection rules discussed above. For a given vibrational transition, Eqs. 3.1/ to 3.19 
describe the line centers of each rotational transition (with Eq. 3.18 not allowed for 
S transitions). Recalling also that the 0-0 rotational transition is always forbidden. 


t/ Bv u ,Av l 


can be evaluated as follows. 
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c 

L 'Bv*,Av l ^\^ Bv u tAv t 2J i max + 1 
where the first form is for non-E transitions, and the second for E transitions. 
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Substituting Eqs. 3.17 to 3.19 in Eq. 3.40 and 3.41 and collecting terms results in 
the following expressions for the averaged frequency. 
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(3.43) 


where again the first form is for non-E transitions, and the second for E transitions. 
While these equations could be implemented as is, it is simpler and still accurate to 
again neglect the D e terms compared to the B v terms. The result is: 
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[3.45] 


These equations for non-E and E transitions can then be substituted in Eq. 3.30. and 
the averaged spectral absorption coefficient can be calculated. 
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The final result for the '‘smeared band ' molecular absorption coefficient is 


= q "“' 
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(3.46) 


where F is given by Eq. 3.44 for non-E transitions, and by Eq. 3.45 for E transitions, 
and the vibrational energy E v t is evaluated from Eq. 3.10 with the appropriate values 
for the lower energy level. To be consistent, the approximation 2J l + 1 » 2 J l is used 
to reduce the fraction inside the summation over J . 


Maximum Rotational Quantum Number 

The maximum rotational quantum number J max can be determined to varying 
levels of approximation. The first order calculation would predict the maximum to 
occur when the energy of a level exceeds the dissociation energy of the species. As 
discussed in Herzberg [65], however, predissociation may occur due to the shape of 
the energy potential curve. Whiting et al. [19] developed a computer code based on 
this fact to predict the maximum rotational quantum number. The method was later 
refined by Whiting [68]. This most recent method has been adopted for use here. 
For each vibrational band of each molecular contributor, it computes the maximum 
allowable rotational quantum number to be considered, based on the shape of the 
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Morse-centrifugal potential describing the rotating and vibrating molecule. 

Further enhancements to the theory for the maximum rotational quantum number 
are possible, but it is felt that this level of accuracy is sufficient for the current work. 

3.2 Emission 

Under equilibrium conditions, the emission coefficient is found from the absorption 
coefficient using Kirchhoff's law. 

j*(s) = Ms)/„ P (r)(l -exp(-hvfkT)) (3.47) 



which depends on the existence of an equilibrium temperature T . 

In nonequilibrium, the single temperature T is not defined. Two approaches for 
finding the emission coefficient are possible. A quasi-equilibrium or excitation tem- 
perature appropriate to each radiative transition can be calculated, from which indi- 
vidual contributions to the emission can be obtained. This is the approach adopted 
in the NEQAIR code [20]. Alternately, expressions for the emission coefficient can 
be found for nonequilibrium from considerations of detailed balancing. In principle. 
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these two approaches should lead to the same result. However, the definition of an 
excitation temperature for use in the first approach is very sensitive to predicted 
nonequilibrium energy level populations that depend on rate data with significant 
uncertainties. An excitation temperature is defined by assuming that a Boltzmann or 
Saha equilibrium exists between the upper and lower level populations for a transition. 
For some nonequilibrium population distributions negative excitation temperatures 
can be computed. For bound-free transitions, for example, solving Saha’s equation 

for T gives: 

/ — E n n 

Tex= k\n{NJN + /NJf(T e )) 

If the nonequilibrium population of level n is small the argument of the logarithm may 
be less than one (the function f(T e ) in this equation is a combination of partition func- 
tions), resulting in a negative excitation temperature. Such negative temperatures 
have no physical meaning and therefore require that some justification be developed 
if they are to be used. The detailed balancing approach avoids this difficulty and is 
therefore employed in LORAN. The necessary expressions for the emission coefficients 
of the various radiative transitions are developed below. 

3.2.1 Atomic Mechanisms 

Bound-Bound Transitions 

The radiative emission from an atomic line is given by. [6] 


J e = \ u hu ul Aui 


( 3 . 51 ) 
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where A u i is the Einstein coefficient for spontaneous emission of the transition from 
the upper atomic energy level u to the lower level /. Detailed balancing relates A u , 
to the absorption oscillator strength f lu as follows: 


. _ 8jtV a 2 

Alu — ”3 ^ulJl u 

mcr g u 


(3.52) 


Combining these two equations leads to the spectral emission coefficient for bound- 
bound atomic transitions. The line shape 6 is the same as that for absorption. 


= 2 J^2LvlJ lu N u b(v,N e J t ,T e ...) 


me? g u 


(3.53) 


Bound-Free Transitions 

Bound-free emission results from the capture of a free electron into an ion. The 
cross section for this process is given by Zel’dovich and Raizer: [6] 

_ 128tt 4 Z 4 e 10 (3.54) 

<Tcn 3^/3 mc?h 4 v 2 n 3 v 

where v is the initial speed of the captured electron. The number n c of such captures 
into the nth energy level for electrons with speeds between v and v + dv per unit 

volume per unit time is 

n c = N+N e f{v)dv v a cn (3.55) 

The electron speeds are distributed according to the distribution function /(v) which 
is assumed to be a Maxwell distribution in equilibrium at the electron translational 
temperature T e . The energy emitted, hv ul , where the subscript u now denotes the 


free state, is given by the sum of the initial kinetic energy of the electron and the net 
energy required to reionize it from the final energy level /. 


mv 2 . , _ v 

hv u i = — p + (/ - Ei) 


( 3 . 56 ) 


I is the ground state ionization potential. Ei is the energy of the bound level, mea- 
sured from the ground atomic state. This equation is solved for v 2 and substituted 
to eliminate v in favor of the frequency v in Eqs. 3.54 and 3.55. 

Combining all these pieces, the emission coefficient for bound-free transitions is 
found from: 


]i<lv = Yh j-hi/ u in c dv 
final levels 1 4x 

The complete expression is 
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(3.58) 


The lower limit on the summation corresponds to capture of a zero energy electron, 
nmin is determined for each frequency by setting v to zero in Eq. 3.56. The upper 
limit is determined by ionization of the atom. 


Free-Free Transitions 


Free-free transitions are a special case in nonequilibrium because they involve only 
free electrons. Since electrons are assumed to equilibrate rapidly to some temperature 
T e , Kirchhoff’s law and the development of the emission coefficient lead to the same 
result. 


8/ 27 r \ l / 2 Z 2 e 6 , . 


( 3 . 59 ) 



3.2.2 Molecular Mechanisms 


A "smeared band” expression can be developed for the molecular emission coeffi- 
cient by a process parallel to that for absorption. The total emission in a rotational 
line, J*,, is obtained from the Einstein coefficient for spontaneous emission: 


J u l — A u d a (/lt' u / 


3.60) 


The emission intensity at a particular frequency is obtained by dividing by the solid 
angle and introducing the line shape function, F{v). 

(3.61) 
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Substituting Eq. 3.20 for the Einstein coefficient, this becomes 


lg^.3 

jt ul = N, 3c2 -^ D] iba q v u v i pjvji F(v ) 


(3.62) 


To obtain the “smeared band” result, this expression must first be summed over J 1 
and the average frequency F Bv u Av i introduced. Using the summation rule for pjuji 


this gives 


10^3 

il . = D 2 elBA q v u v i F(v) 


(3.63) 


Proceeding in parallel to the development for absorption, the sum over v is now 
carried out. 

i £ —3 co 

(3.64) 


16x 3 o . v — ' 4 

jl A , Bv uju = D elBA 


where again the frequency remains inside this summation for accuracy. The average 
emission coefficient for a single rotational line is now obtained using the mean value 
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theorem. 
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where A v is the width of the line. 

This equation can now be solved for the nationally averaged emission coefficient. 
The complete "smeared band” emission coefficient is then obtained by carrying out 
the summations over the final rotational and vibrational levels, J u and t’ u . The result 


is: 
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(3.66) 




As in the case for absorption, the infinite sums are truncated by the occurrence of 
dissociation. 

The average frequency, F Sv u,^, in Eq. 3.66 is the same as that for absorption, 
and is given by Eqs. 3.44 and 3.45 for non-E and E transitions, respectively. The 
average line spacing is also the same, and is given by the expression in Eq. 3.37 with 
the rotational quantum number expressed in terms of the wave number according to 
Eq. 3.38. The final term in the emission coefficient, the nonequilibrium upper le\el 
population, is found from the expression: 


Nbv u j u — Nb 


exp( — E v u/kT u ) gj» exp ( — E^JJcTr) 


(3.6' 


Qv{Tv) Qr(Tr) 

where the vibrational and rotational energies are given as a function of wave number 
by Eqs. 3.10 and 3.39, evaluated for the upper energy level, and the electronic le\el 
population, N b , is obtained from the QSS method. The evaluation of the partition 
functions was discussed in the section on absorption. 
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Substituting all these expressions in the emission coefficient, Eq. 3.66, the "smeared 
band” result is obtained: 
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(3.68) 


Again the approximation 2 J“ + 1 « 2J“ is used inside the summation over J“, to be 
consistent with the level of accuracy of the result. 


3.3 Induced Emission 


The above expressions for the emission coefficient do not include all the radiation 
emitted by the medium. Stimulated or induced emission also occurs due to the 
presence of photons. This emission is proportional to the radiative intensity, /^, and 
is therefore commonly included as a correction to the absorption coefficient. The 
induced emission is given by: 



(3.69) 


In equilibrium the intensity, in this equation becomes the Planck black body 
function, and the result of Eq. 3.47 is obtained. The absorption coefficient corrected 
for induced emission can then be identified in equilibrium from Eq. 3.47 as 


= Ms)U — exp (—hv/fcT)) 


(3.70) 
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In nonequilibrium, the intensity cannot be replaced with the Planck function. In- 
stead, the induced emission coefficient of Eq. 3.69 must be substituted in the transport 
equation, along with the coefficients for spontaneous emission and induced absorption. 
The transport equation, Eq. 2.4, then becomes: 

a/ " ( T- ill + „„(*)/,(>, n, i) = ji(s, t) + j Us, t ) (3.71) 


where j* is the spontaneous emission coefficient. Substituting Eq. 3.69 and gathering 
terms containing /„ then allows the definition of a corrected absorption coefficient 
equal to the factor multiplying /„. The result is: 
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/ 

1 / 


= K u 



(3.72) 


The corrected absorption coefficient may be negative in nonequilibrium. This fact 
must be carefully considered when developing any solution algorithm for radiative 
transport in nonequilibrium problems. 

Negative values of k' u occur when the nonequilibrium populations are such that 
the second term on the right-hand-side of Eq. 3.72 is larger than the first term. 
Physically, this means that the upper energy level of a transition or transitions is 
overpopulated relative to the lower energy level. Such a population inversion can 
occur in a nonequilibrium boundary layer, for instance, when the higher energy levels 
are still populated according to the temperature and chemistry of the inviscid shock 
layer. Because of the v~ z term appearing in Eq. 3.72, negative values of tend to 
occur at the low energy end of the spectrum (see Fig. 5.3). 
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4 Transport Theory 

Solving the equation of radiative transport (Eq. 2.4) in a tractable computation 
requires a number of approximations. It is not possible to consider an infinite number 
of frequencies or an infinite number of directions. Assumptions commonly made are a 
step-wise variation of radiation properties over part or all of the frequency spectrum, 
and a one-dimensional variation of the gas properties. 

The first assumption may be acceptable in many situations, but must be used 
cautiously in nonequilibrium. The absorption coefficient in nonequilibrium cannot 
be found simply from temperature and equilibrium composition but depends on the 
population of excited states, and hence the chemical and excitation kinetics. This 
introduces a number of new variables and makes the design of a step absorption 
coefficient difficult. Some method of selecting appropriate frequencies at which to 
compute the radiative properties must be devised. This problem will be discussed in 

Ch. 5. 

The second assumption of one- dimensional property variation should be nearly 
valid in the stagnation region of flows over blunt bodies. However, the results shown 
in Figs. 6b and 11 of Candler [4] suggest that even near the stagnation point this 
approximation may be in error. The assumption certainly introduces considerable 
error in the radiative flux at corners, and is not at all valid in the wake region of 
vehicles such as the proposed ASTV. Such wake flows are of great interest, because of 
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the desire to place unshielded payloads in the lee of these 


vehicles for orbit transfer. 


4.1 Flowfield Coupling 

Including gas radiation effects in a thermochemical nonequilibrium flowfield cal- 
culation means including the coupling between the fluid dynamics and the radiation. 

The coupling occurs because of the radiative flux term V • & ( denoted Q”* b >' Gn ° ff ° 
[69]) which appears in the total and vibrational-electronic energy conservation equa- 
tions. In the tangent slab approximation, the divergence of the radiative flux reduces 
to simply the normal derivative of <?r. The appropriate term is obtained from the 
radiative transport solution as discussed below for each solution method considered. 


4.2 Plane-Parallel Medium 

A common approximation for solving the radiative heat transfer equation is to 
assume that the medium through which the radiation travels varies in only one direc- 
tion. This is often a reasonable assumption in the stagnation region of a flowfield. as 
well as in other geometries which are of interest for radiant heating. In this situation 
the path variable * can be replaced by the Cartesian coordinate y using the chain 

rule: 

d _ d dx_ 

ds ~ 'dxls + dy ds dzds dy ds 


(4.1) 
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Figure 4.1. Shock Layer Geometry for Plane- Parallel Approximation 
For a shock layer geometry, if both s and y are taken along the surface normal direction 

as shown in Fig. 4.1 then 

= ± (4.2) 

ds dy 

Equation 2.4 then becomes 

+ <(!,)/,(», () = JUt, i) ( 4:3 > 

dy 

This is the equation of radiative transfer along the y-axis in a plane-parallel medium 
under conditions of nonequilibrium. 

4.2.1 Transparent Gas 

A transparent gas is one in which the absorption coefficient is assumed zero. In 
this case, Eq. 4.3 can be further simplified to obtain: 

dU . y - - = jl(y, 0 (4 ' 4) 

dy 

This equation can easily be integrated to find 

I„(y, t) = f V jl(y, t) dy 

Jyo 


(4.5) 


For a solution over a flowfield grid, in which the emission coefficient is known at a 
number of discrete points, this equation can be solved by any numerical quadrature 
method. For a transparent gas, the radiative flux <f R can be found simply by multi- 
plying the intensity by 2 7r. The divergence of <f R , needed for coupling to a flowfield. 
can be obtained by numerical differentiation. Despite the assumption that the gas 
does not absorb radiation, this term will be nonzero due to emission from the gas. 


4.2.2 Absorbing Gas 


For an absorbing gas Eq. 4.3 can be integrated formally to obtain an expression 
for the radiative intensity for any location y. This has been done in equilibrium 
for instance by Ozisik [21, p. 267]. It is convenient to split the intensity into two 
components: one directed along +«/, denoted /+, and the other along —y, denoted 
j— ' Ozisik’s equations can be adapted to the nonequilibrium situation by recognizing 
that the source function, S„, in a non-scattering nonequilibrium gas becomes ;'£/<,• 
The formal solution for these two components of the radiative intensity can then be 
written for the nonequilibrium problem as: 


V ^ e-^-Vdrl 

(4.6) 

k 

(4.7) 


where r„ is the optical variable defined for each frequency v as 
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The = 0 boundary is the vehicle surface, which is assumed to emit radiation 
according to Planck’s function at the wall temperature T w . A radiation equilibrium 
wall temperature can be computed for this purpose, but in this analysis T w is assumed 
to be given. The r„ = r 0 „ boundary is the freestream gas ahead of the vehicle. The 
radiative intensity emanating from the freestream is assumed to be zero. The second 
equation therefore reduces to: 


Nicolet [24] assumed a log-linear variation of properties to solve these equations. 
This formulation does not accommodate negative values of the absorption coef- 
ficient corrected for induced emission, which can occur in nonequilibrium. Nicolet s 
transport algorithm can therefore not be used here. As in the transparent gas case, 
any numerical quadrature method can be used to obtain the radiative intensity dis- 
tribution along the body normal, but it must allow for negative <• 

To obtain the radiative flux in an absorbing, plane parallel medium it is easy to 
show that the above equations apply with the simple replacement of by 2r„ and 
the addition of a factor of r. The net radiative flux is then given by 


, , / ,+ , m -2 tv , f Tu Jt(K) 2(r u -rl) d(2T >) _ P" 

9^) “ * + J 0 K'K) 6 1 J Jr, <(K) 


(4.10) 


The divergence of the radiative flux is obtained by differentiating this equation with 
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respect to y , then integrating over the complete spectral range. The result is: 

| dv 

(4.11) 

where t„ is a function of y. For a solution including all the radiative energy, r'min 1® 
zero and t/ max is infinity. In practice, finite limits on the spectral integration will be 
required. Most gas radiation in shock layers in air occurs at energies below 16.5 eV. 
An alternate upper limit at 6.2 eV includes only the visible and infrared portions of 
the spectrum. This range of energies is easier to measure experimentally and includes 
minimal self-absorption. The lower limit is often set around 0.3 eV because a zero 
energy photon causes a singularity in the free-free absorption coefficient (Eq. 3.9). 
The singularity arises because additional quantum mechanical considerations must 
be included at very low energies. The radiative energy omitted by ignoring these low 
energies is minimal. 

4.3 Multi-Dimensional Medium 

Since an exact solution of the radiative transport equations in multi-dimensional 
media is not feasible, an approximate method is sought. Such a method should intro- 
duce considerable simplifications to the governing equations in keeping with the stated 
objective of providing a relatively fast solution. Numerous approximate methods are 
available in the literature. A brief summary was provided in the topical review in 
Ch. 2. The best-developed of these methods are the moment methods, which reduce 
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the problem from a solution of integro-differential equations to a solution of simpler 
ordinary differential equations. Many of these methods have been developed specifi- 
cally for one-dimensional radiative transport. Sparrow and Cess [70, Sec. 7-9] caution 
that the extension to three dimensions is open to question and must be made with 
caution. They also show, however, that the moment method reduces to the correct 
expressions in the limit of optically thin or optically thick gases and suggest that it 
should therefore be of reasonable accuracy for all values of optical thickness. Modest 
has modified the differential approximation to make it even more accurate for small 
optical depths and generalized it for the three-dimensional case. His method gives 
results of excellent accuracy for all optical conditions [27]. It is adapted below for the 
current non-gray, nonscattering, and nonequilibrium flow. 

As previously mentioned, the source function in a nonequilibrium, nonscattering 
medium is 



The definitions of the incident intensity and flux are repeated here for convenience: 

G,(r)= / I u {f,(l)du (4.13) 

q. R (f) = (4.14) 

In what follows the R subscript on will be dropped. The intensity is now divided 
into tw'o separate contributions: one due to the emission from the medium, / m : and 
the other traceable to the wall emission, I w : I{r,Q) = I w {r,Q) 4- I m (f.Q). The 
incident intensity and flux will then also have two components: G = G w + G m and 
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q = q w + q m . While the intensity attributable to the wall may be highly directional, 
that emanating in the medium should vary quite smoothly. This portion is therefore 
assumed to be given by the P-1 approximation, which is the simplest form of the 
moment method. In this case, the medium intensity can be expressed as 

I m * ^~{Gm + 3<f m • fi) ( 415 ) 

This expression is unchanged from the work of Modest [27]. The equation governing 
the medium intensity then becomes 

s7-(nhJ = jl-< I ^(r, n) < 4 * 16) 

where the nongray, nonscattering, nonequilibrium nature of the gas has been ac- 
counted for. Taking the zeroth and first moments of this equation and integrating 
over 4tt steradians yields the applicable governing equations for G v and <f„ of the 
medium: 



(4.17) 

V 

(4.18) 


The details of the numerical solution to these equations are given in Appendix B. 

To obtain the divergence of the radiative heat flux, which is required to couple 
the solution with the flowfield, the contribution to q from the walls must also be 
considered. This is obtained from the following integral. 

q Va = - I B l/ [r(u>)]e~ Tl/ tldu) 

7T JA x 


(4.19) 
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where the wall radiosity, B u , is given by the sum of emission at the wall temperature 
T w , and reflection of radiation emanating from other wall surfaces. In Modest's 
analysis, the wall reflection is assumed to be diffuse. The treatment of specularly 
reflecting wall surfaces adds considerably to the complexity of the problem and will 
not be addressed here. In any case, for a forebody flowfield the wall surface is convex 
so that no reflection of radiation emitted by the wall can occur. Then the radiosity 
is simply the emission from the wall: 

B„(r) = e v xlvp m {r) (4.20) 

This contribution to the radiative flux must be computed for each grid cell in the 
medium by integration over all the wall surface elements to which it has a line of 
sight. If it is further assumed that the wall is cold, or that it makes a negligible 
contribution, then q w ~ 0, and only the medium contribution need be considered. 
This approximation is not unreasonable for aeroassist flowfields, since heating rates 
are relatively low in the upper atmosphere and the wall temperature will accordingly 
be relatively low during much of the flight. In fact, at the maximum wall temperature 
predicted for the Aeroassist Flight Experiment (AFE), the peak of the black body 
emission spectrum occurs at about 0.75 eV. This temperature is close to the limit 
achievable for reusable surfaces using Space Shuttle-type thermal protection system 
tiles, and so is representative of the maximum wall temperature to be expected for 
a non-ablating surface. For lower temperatures, the peak shifts to even lower ener- 
gies. In this low energy spectral region the gas is quite transparent to radiation (see 
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Fig. B. 1 ) so the wall can have little effect on the divergence of the radiative flux. In 
that case. V • * * V • fa and is obtained by integrating Eq. 4.17 over frequency. 
Note that unlike the plane-parallel medium, numerical differentiation is not requtred 

here. 

Boundary conditions in the multi-dimensional case are less obvious than for the 
plane parallel approximation. It may again be assumed that the radiative intensity 
emitted in the freestream is zero. Even in a three-dimensional flow-field, however, 
symmetries are generally exploited. This means that some grid boundaries will be 
symmetry planes. Others will be outflow boundaries. Symmetry planes can be treated 
with reflecting boundary conditions, but the treatment of outflow boundaries is til- 
defined. The conditions of the gas are not known beyond such boundaries, though tt 
is expected to continue radiating for some finite distance for most flight conditions of 
interest. The treatment of outflow boundaries is discussed further in Sec. 6.2.1 and 

Appendix B. 

For the medium intensity, the vehicle surface is considered to be a cold wall, for 
which Marshak's boundary condition [21] for the P-1 method is 

■2 [1 - l] £. • n + = 0 ( 4 ' 21 > 

The flux, if any, impinging on a wall surface due to emission or reflection 
wall surfaces must be obtained by an integration which considers the view factors of 
the particular problem. When the wall is convex, this contribution is zero. The nu- 
merical details of the treatment of the boundary conditions are given in Appendix B. 
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5 Implementation 

The implementation of the theory for radiation transitions and radiative transport 
developed in the previous two chapters in a practical computer code requires some 
additional work. The details of the process are described in this chapter. 

5.1 Coding 

The present method is implemented in the LORAN (Langley Optimized RAdia- 
tive Nonequilibrium) code, which was developed using two existing codes as starting 
points. The RAD/EQUIL program [33] was followed in developing some of the code 
structure. It also is the source of part of the algorithm used to model atomic line 
radiation. The quasi-steady-state (QSS) excitation portion of the NEQAIR program 
[20] has been minimally adapted for use in LORAN. 

LORAN is divided into four major sections: reading data and setting up proper- 
ties which are constant for a calculation, calculating the nonequilibrium excitation, 
evaluating the radiative properties, and solving the radiative transport problem. 

5.1.1 Initial Setup 

Data Input 

To evaluate the expressions for absorption and emission obtained in Ch. 3 for each 
of the radiating mechanisms, input values are required. Data are also required to 
predict the nonequilibrium populations of the various bound energy levels. Flowfield 
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data are required to define the local gas conditions. Finally, program control inputs 
are desirable to maximize the flexibility of the method without requiring constant 
reprogramming and recompiling. These should allow the user to turn individual 
radiating mechanisms or species on or off to perform detailed studies of the radiation 
spectrum. They should also provide for selection of options in the code (i.e. tangent 
slab vs. 3-D transport). These inputs are obtained from a variety of sources, as 
discussed below. Most of the inputs are species properties and will not change unless 
significantly different new data become available. 

Flowfield Data 

The information required about the gas conditions consists of the number densities 
of each of the radiating species considered and temperatures for the several energy 
modes. The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) 
flowfield solver [71] is the source of the flowfield data used in developing this method. 
It is an 11-species model which provides number densities for all the nitrogen/oxygen 
radiating species of interest. LAURA currently incorporates a two-temperature model 
for thermal nonequilibrium [72]. The rotational and heavy particle translational en- 
ergy modes are assumed to be equilibrated, while the vibrational and electron trans- 
lational energy modes are assumed to be in equilibrium with each other but not 
necessarily with the other modes. This means that the separate T v and T e carried in 
the development of the absorption and emission coefficients are equal, as are T r and 
T t . The distinction between these temperatures was maintained in LORAN, however, 
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so that each can be used if the separate temperatures become available from LAIRA 
or another Computational Fluid Dynamics (CFD) code in the future. Some recent 
work [73] suggests further that distinct vibrational temperatures may exist for each 
molecular species. Should a CFD code ever include such detail, minor modifications 
to LORAN will allow these individual T v values to be used. 

An effort has been made to keep the input of the flowfield properties from LAURA 
generic, so that with minor changes to a single subroutine LORAN can be interfaced 
with a different CFD code. 

Maximum Rotational Quantum Number 

The radiative property data discussed above are read into the program first. As 
part of this process, the maximum rotational quantum numbers are assessed using 
the algorithm developed by Whiting (see Sec. 3.1.2). These values do not change in 
the program. 

Continuum Spectrum 

The spectral location of all radiating transitions is determined from the set of 
energy levels in the atoms and molecules appearing in the gas. This information can 
be used to select a minimum set of spectral points for use in the radiation calculation. 

Consider the bound-free continuum. The absorption coefficient for this mechanism 
is a smooth function of frequency with discontinuities corresponding to the activa- 
tion of additional energy levels. To resolve this spectrum, spectral points must be 
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Figure 5.1. Example of Spectrum Optimization for Bound- Free Continuum 

placed within a very small interval on either side of each jump, as shown in Fig. 5.1. 
However, this interval is much smaller than what is needed for the smooth portions 
of the curve. RAD/EQUIL took advantage of this fact by allowing the user to input 
a few spectral points to capture the jumps. In LORAN, this procedure is automated 
and incorporated in the radiation calculation. For the particular set of atomic species 
considered in a computation, a tailored atomic spectrum is generated to resolve the 
active discontinuities, as the figure illustrates. This simplifies the inputs and elimi- 
nates the chances of omitting a level. It also allows the atomic continuum calculation 
to be performed on the minimum resolving array of spectral points. 

As discussed in Ch. 3, the molecular radiation has been modeled here using the 
“smeared band” approximation. In this approximation, the rotational structure is 
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smoothed into an exponential variation. To resolve the resulting spectrum, spectral 
points must be placed around each vibrational band head. A few points are sufficient 
to capture the smoothed rotational variation within each such band. The resulting 
structure is similar to the bound-free spectrum, though it is much more complicated 
because the vibrational bands are far more numerous and often overlap each other. 
These considerations have been used to develop a routine to choose an optimum 
molecular spectrum for each input spectral range and set of active species. 

For each active species, the vibrational band heads that occur in the spectral 
range under study are computed from Eq. 3.36 with J = 0. An initial spectrum is 
generated by resolving the band head and spacing nv points in the interval between 
the zero rotational level (J = 0) and the maximum rotational level ( J = J m ax)- An 
example for a major vibrational band of iV 2 + is given by the curve labeled "initial 
spectrum” in Fig. 5.2. Note that the clustering of points along this curve results 
from additional weaker bands which underlie this major band. Since the absorption 
coefficient decreases rapidly for a given vibrational band (note the log scale in the 
figure), this spectrum can be reduced. This is done by discarding points that describe 
a vibrational band away from its band head (points required to resolve the jump are 
kept), if they are within a given interval crit of another point. This procedure can 
reduce the spectral array by nearly an order of magnitude. A second reduction can 
be made by discarding any point which has a neighbor within a smaller interval 
crit'2. Applying these steps to the spectrum in Fig. 5.2 results in the optimized 
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Figure 5.2. Example of Spectrum Optimization for a Vibrational Band 

spectrum shown. Applying them to the entire set of vibrational bands, a spectrum 
can be generated that resolves the molecular transitions using less than 2000 points 
and that is optimized for whatever species and spectral range have been requested. 
This calculation is computationally expensive, since it sifts through a large number 
of possible spectral points to produce the optimized molecular spectrum, and it is 
therefore only exercised once per run. The resulting spectrum is stored for all further 
computations. Appropriate selection of the parameters nv , crit, and crit'2 controls 
the amount of detail included in the molecular spectrum and its size. 

Once an optimized spectrum is obtained, the index of the first and last spec- 
tral point for each vibrational band is determined. These indices simplify the logic 
required in computing the molecular band radiation, and improve its efficiency. 
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A few (about 50) points from the optimized molecular spectrum that resolve 
strong vibrational transitions are identified and interleaved into the atomic continuum 
spectrum discussed above. The resulting spectrum of about 150 points will be referred 
to as the modified atomic spectrum. Its use is discussed in Sec. 5.1.4. 

5.1.2 Nonequilibrium Excitation 

To obtain the nonequilibrium electronic level populations for the atoms and molecules 
in the flowfield, data on the rates of excitation and deexcitation of each level are re- 
quired. Sources for this data are varied, but the rates presented are often uncertain 
due to the experimental difficulties involved in their measurement. Park [20] has col- 
lected a complete set of rates for nitrogen and oxygen species which has been adopted 
here. 

The nonequilibrium excitation calculation in the present method was taken from 
Park’s NEQAIR program. It consists of a QSS solution to the set of excitation rate 
equations, in which it is assumed that the rate of change in a level’s population is 
small compared to the rate of transitions into and out of that level [67, Chapter 3]. 
This method is designed for use in situations where the amount of nonequilibrium is 
“not too large”. Its most serious weakness in the present application is that it does 
not account for the history of the flow in areas of high gradients (boundary layer and 
shock). 

The excitation calculation depends on excitation rate data for the various energy 
levels. This type of data is rather difficult to obtain and therefore has a significant 
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level of uncertainty. Park [3] cites this uncertainty as a factor of two, and asserts that 
the impact on the excited state populations is negligible. Further verification of this 
accuracy would be desirable in the future. 

5.1.3 Radiation Properties 

Bound-Bound Mechanisms 

For bound-bound transitions the oscillator strengths fi u , line centers v C l and 
Stark (half) half-widths 7 s are needed. The National Bureau of Standards (NBS) 
has published a large number of line centers and oscillator strengths [66], while Stark 
(half) half- width information is available from Griem [63], among others. In addition, 
the energy levels and degeneracies of the upper and lower levels of the transitions can 
be found in the NBS compilation. 

In the RAD/EQUIL code, Nicolet approximated the so-called high series lines by 
an integral. These lines have lower state quantum numbers of four or greater, and are 
numerous and closely spaced at energies below 3.2 eV. Park treats each of these lines 
individually in NEQAIR. In the present method, in order to reduce the number of 
lines computed and therefore the run time, the high series lines were represented by 
multiplet-averaged values. Data for these averaged lines were obtained from the NBS 
tables, for transitions corresponding to Park’s inputs. This approximation should 
have minimal impact on the results, and yet allows the number of lines treated to be 
reduced by more than half. 

The line shape models for LORAN are taken from the development by Nicolet 
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[23] as implemented in the RAD/EQUIL code [33]. This includes natural, Doppler, 
resonance, and Stark broadening of the atomic lines, as discussed in Ch. 3. Lorentz 
or Doppler line shapes are selected depending which line width is larger. Each line is 
resolved by a small number nppl (generally less than or about 15) of spectral points 
that are distributed starting from the line center (see Fig. 3.1). The location of 
the spectral point farthest from the line center is determined from the width of each 
atomic line. In Nicolet’s implementation the line width used to distribute these points 
was an approximation to the maximum width in the gas layer so that all or most of 
the energy in each line would be included. This approximate maximum width was 
calculated at a location input by the user at which the maximum temperature was 
expected. 

For an application in which a complete flowfield is to be computed, the identifi- 
cation of the occurrence of maximum line width presents some difficulties. Further, 
since the line width may vary by orders of magnitude between various points in a 
flowfield, the use of a spectrum geared to the largest width can result in significant 
loss of accuracy. For these reasons, the routine which sets up the line frequency 
spectrum (subroutine freq of RAD/EQUIL) has been modified to be applied to every 
point in the flowfield. This results in the specification of a local set of spectral points 
based on the local line width. These disparate frequency spectra are then reconciled 
to generate a single atomic line spectrum. This process results in the selection of 
a set of spectral points capturing the energy in a line at its broadest, but involving 
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detailed spectral information at each grid point. This technique allows satisfactory 
resolution of the atomic lines, without requiring an excessive amount of storage or 
computation. 

An additional change made to Nicolet’s method arises because Stark broadening 
is not always the dominant effect on line width. In the original RAD/EQLIL. the 
point distribution logic in subroutine freq considered only the Stark line width. In 
the present model, the largest of the Doppler and Lorentz (Stark effect) line widths 
is used. While this is a small change, it should result in more complete coverage of 

the atomic line radiation in the present model. 

The RAD/EQUIL method has also been modified to use nonequilibrium excited 
state populations and electron temperature, where appropriate. The Doppler line 
width remains a function of heavy particle temperature, since it arises from the ther- 
mal motion of the atoms themselves. 

Bound-Free Mechanisms 

Computing the radiation from atomic mechanisms requires knowledge of the en- 
ergy levels of each atomic species. The large number of electronic energy levels in 
Eq. 3.8 must be modeled by a manageable set of inputs. To reduce book-keeping 
and complexity, a common practice is to group neighboring energy levels into one. 
with an appropriate overall degeneracy. This practice has been adopted here. The 
grouped levels from Park’s QSS method [20] are used for nitrogen and oxygen atoms 
(twenty-two levels for atomic nitrogen, nineteen levels for atomic oxygen). The re- 
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quired values for the principal quantum number n and the energy level E n for each 
of these levels is input. 

The bound-free photoionization edge structure is produced by the activation of 
each level in the summation term of Eq. 3.8. This requires determining the lowest 
accessible bound energy level, n*, for each spectral point in the atomic continuum 
spectrum. The ground state cross sections are invoked as constants to avoid the 
errors arising from the hydrogenic model, as discussed in Ch. 3. Each of the atomic 
energy levels is treated individually in the bound-free calculation, but since these are 
grouped levels many more energy levels are in fact approximately included. 


Free-Free Mechanisms 

LORAN considers only the radiation produced by an electron slowed by the pres- 
ence of an atomic ion. The free-free radiation produced by the proximity of a neutral 
atom or a molecule is considered negligible. The free-free contribution is computed on 
the same spectrum as the bound-free radiation. The resolution thus obtained is more 
than sufficient, and allows the two atomic continuum processes to be easily combined. 


Molecular Mechanisms 


To predict the molecular band radiation, spectroscopic constants describing the 
distribution of rotational and vibrational energy levels according to Eqs. 3.11 and 3.10 
must be provided. These can be found in a number of sources, such as Herzberg [65]. 
Bond et al. [74], and many others. These sources also provide the necessary informa- 
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tion on the electronic energy levels and energy level degeneracies in molecules. The 
intensity of radiation depends on the Frank-Condon factors q v u v t and the electionic 
transition moments ZW Data on these quantities may be found in many often 
contradictory sources. The complete set of molecular inputs collected by Park foi his 
NEQAIR code [20] has been adopted here for simplicity. 

The molecular band radiation is computed using the optimized molecular spec- 
trum whose determination was discussed in Sec. 5.1.1. The expressions for the emis- 
sion and absorption coefficients developed in Ch. 3 have many terms in common. 
Advantage is taken of this fact to reduce the computation time required. Since the 
molecular calculation remains one of the more expensive parts of the radiation calcu- 
lation, care was taken in programming this particular mechanism to ensure that the 
advantages of vectorization are exploited to the extent possible. 

5.1.4 Radiative Transport 

Two distinct methods to compute radiative transport were developed in Ch. 4. 
with and without the tangent slab approximation. Both employ an approximate 
treatment of atomic lines based on the method of RAD/EQUIL. A single mean value 
of the continuum (atomic continuum plus molecular band) absorption and emission 
coefficient is computed for each line and added to the detailed absorption and emission 
coefficients calculated at the nppl points describing the line. These mean values are 
obtained by interpolating the information in the optimized continuum spectra, thus 
avoiding an expensive computation at each of the nppl spectral points describing 
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each line. Since each atomic line covers only a very small part of the spectrum 
(see for instance Fig. B.l), this approximation should introduce minimal inaccuracy. 
Radiative transport is then computed at all nppl spectral points describing each line. 
This detailed result can be integrated over frequency to provide an average flux for 
each line if a concise presentation is desired, or it can be presented in full detail. 

The molecular spectrum can be averaged to obtain values at each point in the 
much smaller modified atomic continuum spectrum, so that the number of spectral 
points for which transport must be calculated is significantly reduced. Since these 
approximate results are based on a detailed molecular spectrum of both emission and 
absorption, the prediction of molecular radiation transport should be quite accurate. 
This averaging technique has been used for the full flowfield solutions reported in 
Ch. 6 to reduce the computation time. For those cases in which only the stagnation 
line is computed, the complete optimized molecular spectrum is used in the transport 
calculation. 

Tangent Slab 

The numerical integration of the transport equation simplified for the case of one- 
dimensional radiation (Eq. 4.3) is straightforward, with one exception. The method 
used by Nicolet in the RAD/EQUIL code assumes a log-linear variation of the absorp- 
tion coefficient between grid points [24]. As mentioned in Sec. 3.3, under nonequi- 
librium conditions there is no longer any assurance that the absorption coefficient 
corrected for induced emission, <, will be positive. The log-linear variation assumed 


by Nicolet cannot handle a negative <. For this work, a simple linear variation 
between grid points is used instead. 

Multi-Dimensional Transport 

The details of the modified differential approximation used to solve the radiative 
transport equation for multi-dimensional media are discussed in App. B, along with 
related numerical considerations. When k' u is not positive, numerical instabilities have 
been observed in this method. Further work is needed to determine whether these 
instabilities can be removed. For now, spectral points at which these instabilities 
occur are ignored. Figure 5.3 is a plot of the number of occurrences of a negati\e 
in a test flowfield at each point in the radiation spectrum. It shows that significant 
numbers of negative only occur for a few points in the low energy end of the 
spectrum. The effect on predictions of wall radiative flux from ignoring these few 
spectral points will be quite small (see Ch. 6). 

It should also be pointed out that the occurrence of these negative values of k' u 
is controlled by the nonequilibrium excitation calculation. As mentioned before, this 
calculation contains significant uncertainties. An example of a population distribution 
at a point with negative k' v is shown in Fig. 5.4 (an equilibrium distribution would 
appear as a straight line in this figure). This population distribution is quite odd, with 
neighboring energy levels having populations different by many orders of magnitude, 
and may be due to uncertainties in excitation rates. It is possible, therefore, that 
at least some of the negative k' u values result from uncertainties in the excitation 
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Figure 5.3. Occurrence of Negative k' u in a Test Flowfield 

calculation. Ignoring these spectral points may therefore introduce less error than 
including them. 


b.2 Computational Optimization 

Obtaining radiation predictions for complete flowfields, or coupling radiation and 
flowfield solutions, requires that the radiation calculations be carried out many times. 
In order to make this practical and to minimize the cost, it is desirable to optimize 
the radiation program, and to determine the minimum set of calculations which will 
provide an accurate radiation prediction. Methods of achieving this are discussed 


below. 


67 



Figure 5.4. Nonequilibrium Population of Excited States of 0 for Negative k' u 

5.2.1 Radiation Calculation 

Parametric studies were conducted to determine the spectral resolution required 
for acceptable accuracy in the predicted radiation. Variables whose effect was studied 
include na, the number of points allowed in the atomic continuum spectrum; nppi 
the number of points resolving an atomic line; and the several variables which control 
the generation of the spectrum for molecular band radiation {nv, crit , crit2). 

Tables 5. 1-5.3 summarize the findings of this optimization process. The baseline 
against which the flux calculation is judged (first line of each table) is a detailed spec- 
trum case with na =1500, nppl =15, nv =10, crit =0.00825 eV, and crit‘2 =0.0011 eV . 
While uncertainties in the radiation calculation mean that this prediction is not nec- 
essarily right, such a detailed spectrum does provide the best result possible. The run 
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times are judged against a case with na =100 and all other values unchanged (second 
line of each table). The transport time appearing in column 4 is the time required 
for tangent slab transport. Relative time savings for the modified differential approx- 
imation should be similar. Alternate values tested for these parameters are given in 
the first column of the tables. The designations Cl and C2 refer to combinations of 
the alternate values as follows: Cl has na =75, nppl =11, nv =8, crit =0.0165 eV. 
and crit2 unchanged; C2 is the same as Cl except na =100. Three flight conditions 
identified as Cases A, B, and C, were studied to obtain results over a wide range of 
nonequilibrium conditions. These cases are described in further detail in Ch. 6 and 
7. 

As shown in Tables 5.1-5. 3, the prediction of the radiative flux is affected by less 
than 5 percent for each set of parameters, except for two occurrences in Case C (see 
column 2). These are both attributable to setting na =75. The potential savings in 
computational time are considerable. Based on these results, it is determined that 
the parameter set denoted C2 [na =100, nppl =11, nv =8, crit =0.0165 eV, and 
crit'2 =0.0011 eV) is appropriate for use in the radiation model. Increasing crit 2 in 
line with the changes in the other parameters was found to have no effect, so results 
are not presented here. 

5.2.2 Excitation Calculation 

The last column in Tables 5. 1-5.3 shows the fraction of CPU time used by the 
radiation and transport calculations. The remainder is required to obtain the nonequi- 



Table 5.1. Computational Optimization - FIRE II at 1631 sec 


Case 

Flux 

CPU time 
for Radiation 

CPU time 
for Transport 

CPU 

Fraction 

Baseline 

0.508 

6.116 

1.735 

0.50 

na =100 

+0.8% 

3.561 

0.895 

0.37 

nppl =11 

+0.9% 

-13.5% 

-29.7% 

0.33 

na =75 

+0.8% 

-1.6% 

-2.1% 

0.36 

nv =8 

+ 1.4% 

-0.3% 

-0.8% 

0.37 

crit * 2 

+0.8% 

-0.9% 

-0.7% 

0.37 

Cl 

+ 1.4% 

-14.9% 

-30.8% 

0.32 

C2 

+ 1.4% 

-26.6% 

-29.5% 

0.30 


Table 5.2. Computational Optimization - FIRE II at 1634 sec 


Case 

Flux 

CPU time 
for Radiation 

CPU time 
for Transport 

CPU 

Fraction 

Baseline 

22.50 

10.843 

3.059 

0.69 

na =100 

+2.1% 

6.326 

1.714 

0.57 

nppl =11 

+3.6% 

-13.0% 

-29.3% 

0.52 

na =75 

+2.5% 

-1.1% 

-1.6% 

0.56 

nv =8 

+2.6% 

+0.1% 

-0.1% 

0.57 

crit *2 

+2.1% 

-0.6% 

-0.8% 

0.57 

Cl 

+4.5% 

-14.8% 

-30.7% 

0.52 

C2 

+4.1% 

-26.0% 

-29.3% 

0.49 


Table 5.3. Computational Optimization - FIRE II at 1637.5 sec 


Case 

Flux 

CPU time 
for Radiation 

CPU time 
for Transport 

CPU 

Fraction 

Baseline 

292.9 

10.061 

2.913 

0.73 

na =100 

+ 1.5% 

5.754 

1.617 

0.60 

nppl =11 

+1.9% 

-12.1% 

-28.9% 

0.56 

na =75 

+6.1% 

-1.2% 

-1.1% 

0.60 

nv =8 

+ 1.2% 

-0.1% 

+0.1% 

0.60 

crit * 2 

+ 1.6% 

-0.7% 

+0.2% 

0.60 

Cl 

+6.3% 

-14.1% 

-30.2% 

0.55 

C2 

+ 1.7% 

-24.5% 

-28.8% 

0.53 
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librium excited state populations, and is a significant fraction of the total. The QSS 
algorithm used to perform this calculation was obtained from NEQAIR and is used 
as a black box. While important savings may be possible in this algorithm, no at- 
tempt has been made to achieve them. It is clear that this is a limiting factor in 
the efficiency of the radiation prediction. It is interesting, but not surprising, to note 
that the excitation calculation takes longer for conditions which are further out of 
equilibrium. 


5.2.3 Radiation Subgrid 


Calculating radiative properties at every point on a fiowfield grid is often wasteful. 
The absorption and emission coefficients only change significantly in regions where the 
temperatures or species concentrations have large gradients. Bolz [30] developed an 
algorithm to automatically select a subset of grid points for the radiation calculation. 
He defines a weighting function which can be adapted to the nonequilibrium situation 

as 
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The ' in these equations denotes the partial derivative with respect to 77 , the normal 
to the wall. The average derivative of the species number densities is 

X = t\Kj\ 

j=i 


(5.2) 
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Figure 5.5. Example of Radiation Subgrid Selection 

for s chemical species. Radiation is then computed for kn < k grid points at equal 
intervals of Zk . An example of such a grid is given in Fig. 5.5. Though weighting 
factors are included for each term in Z k there is no obvious reason to emphasize the 
gradient in one variable more than that of another. These factors are therefore all 
set equal to one. This algorithm has been implemented in LORAN as an option. It 
can be invoked when reducing the computational time is important. If a separate 
rotational temperature is available, or if individual species vibrational temperatures 
are computed, additional terms can be added to the function Zk. 

The subgrid algorithm has been assessed using the three nonequilibrium test cases 
of Ch. 7. The effects are summarized in Table 5.4. The subgrid results (columns 3 
and 4) were obtained using the optimized C2 set of parameters discussed in Sec. 5.2. 
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Table 5.4. Effect of Radiation Subgrid Algorithm 


Case 

Detail spectrum 
All grid points 

20 grid points 

15 grid points 

A 

0.508 

-1.4% 

4-6.4% 

B 

22.50 

4-10.7% 

+ 14.2% 

C 

292.9 

4-3.9% 

+2.9% 


while the detailed spectrum result (column 2) uses the Baseline parameter set to 
establish a reference heat flux (this is again a grid refinement test rather than an 
accuracy check, since even the detailed spectrum result relies on uncertain radiation 
and excitation properties). The CPU time required to compute the excitation and 
radiation for a given case is a nearly linear function of the number of grid points, so 
the time savings can easily be assessed. 

For each of the three cases using twenty grid points (out of 64 in the flowfield 
grid) selected with the modified Bolz algorithm provides radiative heat flux predic- 
tions within about ten percent of the reference result, even including the spectral 
approximations introduced by the use of the C2 data set. A subgrid of this size is 
therefore considered adequate, confirming Bolz’ conclusions. 

In an axisymmetric or 3-D shock layer flowfield this algorithm can be applied to 
each normal ray of grid points. This results in a ‘‘patchwork grid, in which cell 
shapes are more irregular than those in a standard LAURA grid. In an extreme case, 
such a grid may destabilize the transport solution. In most cases, however, optimizing 
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Figure 5.6. Patchwork Grid Resulting from Application of Subgrid Algorithm 

Along Each Normal Line 

the subgrid along each normal ray may be acceptable and even beneficial. Figure 5.6 
presents an example of such a grid for a test flowfield. The irregularities appearing 
in the mesh are in fact quite minor. 

To further reduce the required computational time, every other normal line of the 
flowfield grid might be skipped in the radiation calculation. The required properties 
can easily be interpolated for the skipped normal lines. 

5.3 Flowfield Coupling 

The divergence of the radiative flux, V • <?R> which appears in the energy conser- 
vation equations to couple radiation and fluid dynamic phenomena, is computed for 
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each point in the flowfield by LORAN, either directly or by interpolation of results 
on a radiation subgrid. Coupling can be accomplished by alternatively computing 
the flowfield with LAURA and the radiative flux term with LORAN. Because the 
computation of the radiation requires up to three orders of magnitude more compu- 
tational time than one iteration of the fluid dynamics computation, however, it is 
advantageous to utilize a multitasking strategy. Multiple tasks are assigned to the 
radiation calculation while a single task computes the nonequilibrium flowfield. The 
computational domain is partitioned into three to seven subdomains (depending on 
the number of processors available), each assigned to a processor for calculation of 
the radiation field. When run asynchronously with the single flowfield task, the mul- 
tiple radiation tasks allow the flowfield equations access to the latest radiative flux 
terms in shared memory (multiple LAURA iterations still occur during a single radi- 
ation update). This asynchronous strategy is an efficient use of a multiple processor 
machine since no processor remains idle waiting for another task. Preliminary com- 
putational experiments indicate no instabilities in this approach. The effectiveness of 
asynchronous multitasking schemes using LAURA is further discussed in Ref. [75]. 

As a practical matter, it is often more efficient to first obtain a converged flowfield 
solution using LAURA, then ‘‘turn on” radiation and converge the coupled flowfield. 
If the radiation is only a small perturbation to the flowfield this second convergence 
may require only a few additional LAURA iterations. If there are strong radiation 
effects then converging the perturbed flowfield will require many additional iterations. 
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Though the necessary code to couple the LAURA and LORAN solutions has 
been generated, only a very few preliminary results have so far been obtained with 
Bowfield coupling. The cases run to date confirm that this coupling can be done, but 
no significant effects have yet been observed for the cases so far examined. Coupled 


results are therefore not included in Ch. 6. 
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6 Results and Discussion 

Gas radiation is a complex phenomenon, so that true verification of a radiation 
model requires checking it against a large number of flight and ground-based ex- 
periments covering a wide range of conditions. Only a few such measurements are 
currently available. Representative results for these cases are shown here as an initial 
test of the method. All nonequilibrium flowfields used in the verification were gener- 
ated by the LAURA code of Gnoffo [71], in order to eliminate differences caused by 
flowfield modeling. As discussed in the following chapter, the selection of models in 
a flowfield code can have a large impact on the predicted radiation. 

6.1 Comparison to Experiment 

6.1.1 Ground-Based Data: AVCO Shock Tube 

One of the classic radiation measurements is that in the AVCO shock tube [43]. 
Using the gas conditions suggested by Park [72], a predicted spectrum for the emission 
intensity has been obtained for the peak radiation point, and is shown in Fig. 6.1. This 
spectrum can be compared to Figure 5 of the just-cited reference, which is reproduced 
here in Fig. 6.2 and shows the measured spectrum compared to a prediction obtained 
by NEQAIR. The agreement between the measured data and the present result is no 
worse than that shown by NEQAIR and may actually be a little better. The atomic 
lines appearing in Fig. 6.1 are multiplet-averaged, and so should exhibit only a gross 


agreement with the NEQAIR prediction of Fig. 6.2, as is observed. 

It is clear from Fig. 6.2 that the measurements obtained in the AVCO experiments 
had a lot of scatter. It is difficult, therefore, to do more than note the general 
agreement between the exper.mental results and the predictton. Both predictions do 
have significant differences from the AVCO data in a few spectral regions (particularly 
below 0.3 pm). These may be attributable either to contamination of the shock tube 

or to inexact specification of the gas conditions. 

6.1.2 Flight Data: Project FIRE 

The FIRE flight project of the mid-1960’s consisted of two flights simulating reen- 
try of an Apollo type vehicle. The first flight, FIRE I, experienced telemetry problems 
that made data reduction and analysis difficult as well as control problems during the 
second half of the entry that resulted in substandard data. These problems were 
corrected before the FIRE II flight, which provided good measurements of the total 
and spectral radiation at the stagnation point during the forty second entry period 

[76, 77], 

The FIRE flight vehicles had an Apollo-like geometry with a layered beryllium 
heatshield (the second layer had a geometry identical to Apollo). Each of the three 
layers was used up to a temperature limit and jettisoned, resulting in three periods 
of prime data during the entry. The firs, heatshield had a nose radius of 0.935 m and 
a diameter of 0.672 m. Only data obtained during the use of this first heatshield will 
be examined here, since it coincided with the nonequilibrium portion of the flight. 
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Figure 6.1. Calculated Spectrum for the Peak Radiation Point of the AVCO 

Experiment 
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Figure 6.2. Measured Spectrum and NEQAIR Prediction for the Peak Radiation 

Point of the AVCO Experiment 
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The trajectory of the FIRE II flight is reported in Lewis and Scallion [62]. Fore- 
body temperature histories are given in Cornette [78]. The instrumentation used for 
the flight measurements is described by Richardson [61]. In particular, the spectral 
response of the radiometer windows is reported by Dingeldein [79, Fig. 5]. It is flat 
between 0.23 and 2. p in wavelength, falling off sharply below 0.2 fi m and more 
gradually above 2. fim. The design goal is reported as a flat response from 0.2 to 
6. fim. The various reports on FIRE quote actual spectral ranges starting between 
0.2 and .23 p, and ending between 4. and 6. fim. Taking these variations into 
account, the spectral range of 0.23 to 4. fim (0.31 to 5.4 eV ) is considered to be 
closest to the actual window transmission range. FIRE also had scanning spectral 
radiometers that were designed to cover the range of 0.2 to 0.6 fim (2.1 to 4.1 eV ). 
Mechanical problems during the flight limited forward scans to 0.3-0.5575 p, and 
backward scans to 0.6090-0.3 fim [77]. The spectral resolution is quoted as 0.004 fim , 
with a root-sum-square uncertainty in the measured spectra of ± 23 percent. This 
resolution is insufficient to resolve any atomic lines. 

For this study of stagnation line radiation, the FIRE vehicle geometry has been 
modeled by a sphere with an effective nose radius of 0.747 m. A fully catalytic 
wall boundary condition at the flight-measured wall temperature was used. The 
cases selected from the FIRE II flight are shown in Table 6.1 [62]. The designations 
Case A, B and C will be used below to identify each case. These cases cover the 
range from extreme nonequilibrium to near-equilibrium according to Cauchon [76. 
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Table 6.1. Selected FIRE II Flight Conditions 


Case 

Time 

(sec) 

Altitude 

(km) 

Velocity 

(km/sec) 

Density 

(kg/m 3 ) 

T 

-i- XU 

(K) 

1 

9R 

(W/cm 2 ) 

A 

1631 

84.59 

11.37 

9.15e-6 

460 

.16±20% 

B 

1634 

76.41 

11.36 

3.72e-5 

615 

8.2±20% 

C 

1637.5 

67.04 

11.25 

1.47e-4 

1030 

81.7±20% 


p.H). Predicted temperature profiles along the stagnation line, using the baseline set 
of energy exchange models in LAURA (defined in Ch. 7), are presented in Figs. 6.3- 
6.5. For Case A, extreme thermal nonequilibrium is predicted, in agreement with 


Cauchon, while Case C is predicted to be in thermal equilibrium through about half 
the shock layer. 

The flight radiation levels (,„) quoted in the above table were taken from Cauchon 
[76, fig. 13] and converted from intensities to fluxes assuming a transparent plane- 
parallel shock to be consistent with the calculations below. Also shown is the ro 
sum-square uncertainty estimated in a post-flight error analysis. It should be noted 
that the uncertainty in the radiation level for Case A is probably higher, as the 
radiometers were at the lower limit of their sensitivity range for that case. The 


measurements show considerable scatter during the early part of the trajectory [761, 
with a variation of about a factor of three between high and low readings for the 
first several seconds. The scatter decreases as the level of radiation increases. It is 
noteworthy that the FIRE I data acquired in this density range for a slightly highei 
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Figure 6.5. Predicted FIRE II Stagnation Line Temperature Profiles ■ 163 1 . 5 sec 

velocity (11.57 vs 11.37 km/sec) were higher by a factor of three than the FIRE II 
results for both the total and spectral radiometers. The two flights thus suggest a 
very sensitive dependence of radiative heating on freestream velocity under these flow 
conditions. At the conditions of Case B. the flight-to-flight variation is about a factor 
of 1.6. For Case C, it is about 1.3. 

Detailed comparisons between the LORAN predictions and the FIRE flight ra- 
diometer measurements are given in the following chapter, as part of a study in which 
the impact of various nonequilibrium flowfield models on radiation is assessed. Some 
additional comparisons between LORAN and other nonequilibrium radiation predic- 

tion methods are given here. 

The predicted profiles of radiative emission along the stagnation line for Case B, for 
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Figure 6.6. Emission Profiles for FIRE II at 1634 sec with Comparison to 

NEQAIR 


the spectral range in which the radiometer windows transmit, are shown in Fig. 6.6. 
Calculations from LORAN and NEQAIR [20] are compared, broken out into the 

radiating mechanisms. The NEQAIR result shown here was computed using . 


various 


a version of the NEQAIR code which was obtained from Chul Park in 1988. It has 
been slightly modified to correct some minor errors this version contained, and to 
allow emission calculations within a limited spectral range. The agreement between 
the two predictions is excellent except for a discrepancy in the atomic line radiation 
near the wall. The predicted wall radiative heat flux differs by only 3 percent between 

the two codes. 

Figure 6.7 shows the emission profiles for Case C. The overall agreement be 
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Figure 6.7. Emission Profiles for FIRE II at 1637.5 sec with Comparison to 

NEQAIR 

the two predictions is quite good, although there are some differences. The radiative 
flux reaching the wall varies by about 6 percent between these two predictions. 

Radiation predictions for these cases have also been made by Carlson and Gaily 
[35] for a nitrogen freestream. Figure 6.8 can be qualitatively compared to Fig. 6.9. 
which is reprinted from Figure 5 of their paper. It shows the approximate spectral 
variation of the wall radiative flux over the entire spectrum including the ultraviolet. 
Though the predicted magnitudes are lower in Ref. [35] because only nitrogen is 
considered, the variations are nonetheless qualitatively very similar. The details of 
the spectra cannot be compared because of the line group presentation used in the 
reference. It should also be noted that the chemical kinetics and energy exchange 



85 


Table 6.2. Radiative Heat Flux Predictions for FIRE II 


Case 

NEQAIR 

LORAN 

Difference 

LORAN 


Flux, W/cm 2 (0.31-16.5 eV) 

Percent 

with Absorption 

1634 sec 

1682. 

1806. 

7.4 

29.35 

1637.5 sec 

13233. 

14771. 

11.6 

421.2 


models used by Carlson and Gaily are different from those in the LAURA code. As 
demonstrated in Ch. 7, these models can have a large impact on the magnitude and 
spectral distribution of the predicted radiative heating. 

Results have also been computed for the complete radiation spectrum (0.31 to 
16.5 eV) for these two cases. These are summarized in Table 6.2, which compares 
the emission predictions between NEQAIR and LORAN. The difference, about ten 
percent, is minor. The last column shows the total radiative flux reaching the wall 
when absorption is included in the transport calculation. The large self- absorption of 
the ultraviolet portion of the spectrum is evident here. This wall flux is still somewhat 
high relative to the flight data, but the effect of radiation cooling on the flowfield has 
not been accounted for in these calculations. 

6.2 Nonequilibrium Test Cases 

6.2.1 Mars Return 

A flowfield solution has been obtained for one of a number of possible flight condi- 
tions identified for the return from a mission to Mars [80]. It consists of a 60° sphere 
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Figure 6.8. Spectral Variation of Stagnation Point Radiative Heat Transfer for 

FIRE II 



Figure 6.9. Spectral Variation of Stagnation Point Nitrogen Radiation for FIRE II 
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cone with a 1.08 m nose radius flying at 80 km altitude with a velocity of 12 km/sec. 
The radiative heating for this case is no. too severe, so it provides a good initial 
test of the method. Results have been obtained for this case using both tangent slab 
transport and the modified different, al approximation (MDA). Figure 6.10 compares 
the wall radiation flux predictions from the two transport methods. As expected, the 


MDA result is lower than the tangent slab result. The low *. value for the MDA 
result near the stagnation point is believed to be erroneous, however. The sudden 
increase in at the shoulder in the MDA result certainly is, and is attributable to 


the treatment of the outflow boundary. Further work will be required to correct both 
of these problems. The qualitative trends of the two methods are similar, showing 
the radiative flux increasing with distance from the stagnation point. This increase 
is unlike the behavior of convective heating on such a body and results because ra 
diation is a volumetric quantity. As seen in Fig. 6.12 below, the standoff distance is 
much larger on the flank than at the stagnation point, and the region of significantly 
radiating gas has increased greatly. The vibrational temperature in this nonequilib- 
rium flowfield is still high in this region, and a large amount of radiation is therefore 

emitted from the gas. 


The other variable of interest for coupled flowfields is the divergence of the ra- 
diative flux in the shock layer. Figures 6.11 and 6.12 show contour plots of V • ft 
computed with each transport method. The tangent slab result shown in Fig. 6.11 is 
in fact just the derivative of along each normal grid line, so the values at neigh- 
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Figure 6.10. Wall Radiative Flux for Mars Return Case 


boring grid lines do not influence each other. Compared to the MDA result given 
in Fig. 6.12 the tangent slab result is very disjointed. Accounting for the long range 
influence of radiation with the MDA solution provides a much smoother variation of 
V • ?r. Should stability be an issue in a coupled flowfleld solution, therefore, the MDA 
result might be expected to have less of a destabilizing effect than the tangent slab 
result. While this has not traditionally been used as an argument against the tangent 
slab approximation, it may prove to be decisive at some flow conditions. For the few 
coupled cases run to date, the stability of the CFD solution has not been reduced by 

adding radiation. 


e 6.11. Radiative Flux Divergence 


Case with Tangent Slab 


for Mars Return 
Transport 



, re 6.12. Radiative Flux Divergence for Mars Return Case with MDA 


Transport 
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6.2.2 Aeroassist Flight Experiment 

The Aeroassist Flight Experiment (AFE) is a NASA project intended to obtain 
flight data for the design of future Aeroassist Space Transfer Vehicles (ASTY). Mea- 
surement of the nonequilibrium radiative heating to the AFE surface is one of the 
principle objectives of the project. The aerobrake configuration is a raked cone with 
a blunted elliptic nose [81]. Several studies of nonequilibrium radiative heating for 
this configuration have been made to provide inputs to the design of the heatshield 
[17, 37, 82, 83]. 

The LORAN method has been used to obtain a prediction of the level and distribu- 
tion of radiative heating at the peak heating point in the trajectory. An axisymmetric 
LAURA flowfield solution in which the AFE geometry was modeled by a sphere with 
an equivalent nose radius of 2.16 m [84] was generated [85]. The flowfield models used 
to obtain this solution are the baseline models listed in Sec. 7.1, including the use of 
the Park-87 chemical kinetics. The wall temperature was assumed to be constant at 

1750 K. 

Radiation predictions were obtained for this case, again using both tangent slab 
and MDA transport methods. The tangent slab calculation was completed in about 
45 minutes of actual elapsed time on a Silicon Graphics 4D/320, while the MDA result 
required about 12 hours of elapsed time. The convergence of the MDA solution for this 
case was much slower than that for the Mars return case (in fact this case is not quite 
converged), but even this run time is not unreasonable. The radiative heating rates 
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Figure 6.13. Wall Radiative Flux for AFE 


predicted for the AFE peak heating point are given in Fig. 6.13. These predrct.ons 
are in the same range as earlier computations for the AFE radiative heating [4, 86]. 
The MDA prediction is again lower than the tangent slab result, with problems still 
appearing at the grid boundaries. It is expected that running the MDA transport 
to complete convergence for this case would result in closer agreement between the 
two, as would a correction for the depressed MDA value near the stagnation line. It 
should also be noted that only the nose region of the equivalent sphere flowfield has 
been used for the AFE case. Away from the nose the geometry deviates quickly from 
the AFE configuration so the results are not of interest. 

The radiative flux divergence for the AFE flowfield is presented in Figs. 6.14 and 
6.15 for the tangent slab and MDA transport, respectively. Again the smoother 
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behavior of the MDA solution can be noted. As mentioned earlier, the tangent slab 
method employs a numerical differentiation to obtain dq R /dn. This differentiation 
may be adding to the error already incurred by ignoring the variation of the radiation 

properties in the direction parallel to the surface. 

The amount of heating for the AFE vehicle resulting from radiation in the ul- 
traviolet (UV) portion of the spectrum has been the subject of some debate. The 
distribution of radiation predicted by the LORAN method is shown in Fig. 6.16. 
Curves are shown using the tangent slab and MDA transport models. The relative 
importance of different spectral regions varies only slightly between the two transport 
methods. Both curves indicate that just over half of the radiative heating experienced 
by the AFE spacecraft results from the UV portion of the spectrum. This is a signif- 
icant fraction of the total. Radiation in the UV spectral range is highly self-absorbed 
(see for example Table 6.2). If the amount of self- absorption is mispredicted only 
slightly, the wall radiative flux may be significantly increased for AFE. 
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7 Flowfield Model Studies 


Flovvfield studies were performed using the Langley Aerothermodynamic Lpvvind 
Relaxation Algorithm (LAURA) code developed by Gnoffo [71] to investigate the 
sensitivity of the radiation predicted by LORAN. LAURA, like all nonequilibrium 
Computational Fluid Dynamics (CFD) codes, contains a number of semi-empirical 
models governing the exchange of energy between energy modes. LAURA presently 
incorporates the two-temperature model of Park, in which the heavy particle trans- 
lational and rotational energy modes are assumed to equilibrate to a temperature T t 
(the combination is referred to as the translational mode), and the vibrational, elec- 
tronic, and free electron translational energy modes are assumed to reach a separate 
equilibrium with a temperature T v (the combination is referred to as the vibrational 
mode). Since radiation is strongly influenced by the amount of energy available in 
each mode, a study was conducted to assess the impact that uncertainties in these 
energy exchange models have on radiation predictions. 

The flowfield cases chosen for these studies are taken from the trajectory of the 
FIRE II flight project [61, 62]. Widely varying freestream conditions were selected to 
exercise the exchange models over a range of levels of nonequilibrium. A description of 
this flight experiment was given in Sec. 6.1.2. All the cases run for the flowfield studies 
had 64 grid cells normal to the body in the LAURA solution. This is not sufficient 
to resolve the shock, as illustrated for instance in Fig. 6.5. One grid-resolved solution 
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with 128 grid cells was obtained to quantify the effect of this lack of resolution. 
The radiation prediction obtained with this refined grid was only about 15 percent 
different from that with the coarser grid. This is not enough to change any of the 

conclusions drawn from this study. 

All the results presented in this chapter are for the stagnation line of the FIRE II 
vehicle and use the tangent slab radiative transport method. Only the spectral range 
in which the FIRE radiometer windows were transparent (0.31-5.4 eV) has been 

included. 

7.1 Baseline Models 

A discussion of the energy exchange mechanisms and the alternate models used 
for each is given in the following sections. The set of models considered as the baseline 
consists of the most recent recommendations from Park which have been incorporated 
in the LAURA code. They are, for the vibrational-translational energy exchange cross 

section: 

<7 V = 10 _21 (50,000/T t ) 2 m 2 (‘ D 

for the dissociation temperature: 

T d = T t 7 T v 3 (7 ' 2: 

and for the energy exchange in dissociation: 


AE V = E v 


(7.3) 
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, i.- i ™.-n a« Park-87 f87l is considered part of 

In addition, the set of chemical reasons known as Park 1 

the baseline. 

A run identifier specifies the models used for each solution. The baseline 
denoted „7S. The first two characters, .1, denote the first model for vibrational- 
translational energy exchange, Eq. 7.1. The last two, 7 *, denote the powers m the 
dissociation temperature equation. Only alternate models for the energy exchange 
in dissociation and the chemical reaction set are fiagged in the run rdentifier for 


aanalseness. These will be explained as they are used below. 


7.2 Energy Relaxation 


The modeling of relaxation or equilibration between energy modes is an important 
question in noneqnilibrinm CFD. Generally the process of eqml.brat.on 
pletely understood and must be described semi-empirically. Several formulations have 
been proposed for many of these models. Data to evaluate the different f 


are scarce. 


To quantify the influence which the different proposed formulations of the models 

for energy relaxation have on predictions of radiative heat transfer, a parametric 

, c u models The models and formulations 
study was undertaken on a number of such models. 

considered, and the results obtained, are detailed below. 


98 


7.2.1 Dissociation Temperature 

The dissociation temperature, T d , is the rate-controlling temperature for reactions 
involving the dissociation of a molecular species. Several formulations have been pro- 
posed for this temperature, notably by Park and his co-workers [87, 88]. These models 
are empirical and consist of a geometric weighting of the translational temperature 
T t and the vibrational temperature T v as follows: 

T d = T t m T” (7.4) 

The choice of powers m and n adjusts the weight given to each temperature. The 
recent trend has been to increase the weighting of the translational temperature T t 
because of indications that heavy particle collisions are more important than other 
mechanisms in causing dissociation. Three different models with (m,n) equal to 
(.5, .5), (.7, .3), and (1,0) are considered here. The first set is the original Park two- 
temperature proposal, the second the more recent model, and the third is included 
to study the extreme case of dissociation controlled by the heavy particle tempera- 
ture alone. The dissociation temperature models will be denoted by 55, 73, and 10, 
respectively, in the third and fourth digits of the run identifier. 

As shown in Figs. 6. 3-6. 5, the translational temperature T t is much higher than 
T v near the shock. Increasing the weighting of the dissociation temperature, T d , on 
T t in Eq. 7.4 therefore results in faster dissociation of the molecular species. This 
dissociation removes both translational and vibrational energy, decreasing T) and T v . 
The increased dissociation is shown in Fig. 7.1 which compares the N 2 and 0 2 profiles 
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along the stagnation line for Case A for the two extreme T d models. The effect on the 
temperature profiles is shown in Fig. 7.2 also for Case A. The decreased temperature 
results in increased density and a smaller shock standoff distance. For the three 
FIRE II cases, placing increasing weight on the translational temperature T, in the 
definition of the dissociation temperature T d decreases the total radiation. Par. of 
the decrease occurs in the molecular bands, which depend directly on the molecular 
concentrations. The rest of the decrease results from the lower T u and decreased 
standoff distance affecting all three radiating mechanisms. Figure 7.3 shows the 
emission profiles for the three T d models for Case A. Case A is the most nonequilibrium 
of the three FIRE conditions studied here, and therefore is most sensitive to these 
energy relaxation models. The sensitivity to the T d model decreases at the later, more 

equilibrium, flight times of Cases B and C. 

Tables 7. 1-7.3 at the end of this chapter reveal that in all three cases (vllO, vlSo 
and id 73) the predictions are closer to the flight data when T, is weighted more in 
the definition of T d , confirming the recent work on this model mentioned above. For 
Cases A and B, agreement within about a factor of three is found. For Case C, the case 
nearest equilibrium, the variations caused by this model are within the uncertainty 

of the flight data. 

7.2.2 Vibrational-Translational Energy Exchange 

The equilibration of the vibrational and translational energy modes is modeled 


using a 


relaxation time, r„, for each species 


Millikan and White [89] proposed a 
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Figure 7.1. Effect of T d Models on Molecular Dissociation - 1631 sec 
semi-empirical formulation, t™, in the range of 300 to 8000 K. Park [90] suggests an 
additional collision limiting correction, r p , for high temperatures where the Millikan- 
White correlation predicts excessively fast relaxation. A concise explanation of the 
failure of the Millikan- White correlation is provided by Sharma and Park [91, p. 136]. 
The total relaxation time is then 

r =t mw + t p ( 7 - 5 ) 

• V — * V 1 v 


where Park’s contribution is 

T v = (<7t /njAfj) -1 


(7.6) 


In this expression ~ s is the average speed and ;V 4 the number density of molecular 
species 3 , and is an effective cross section for vibrational relaxation. The latter 
quantity has been the subject of some debate, with three different values receiving 
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Figure 7.2. Effect of T d Models on Temperature Profiles - 1631 sec 



Figure 7.3. Effect of T d Models on Radiative Emission Profiles - 1631 sec 
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support in various efforts to match the limited experimental data that has bearing 
on the question: 


= nr 21 (50,ooo/r t ) 2 m 2 

(7.7) 

<r u = lO" 20 m 2 

(7.8) 

<f v = 10- 21 rn 2 

(7.9) 


Equation 7.8 is the original proposal, which assumes 


this cross section to be one tenth 


of the elastic cross section. Equation 7.9 was suggested when Eq. 7.8 seemed still too 
high. Equation 7.7 is a modification of Eq. 7.9 that reduces the contribution of rf 
at low temperatures. All three models are considered here, and are identified by the 
notation vl, » 2 , and vS, respectively, in the first two digits of the run identifier. 

A further modification to this energy exchange process has been proposed by- 
Park [87], He argues that vibrational relaxation exhibits a diffusion-like behavior at 


high temperatures which requires a 


correction to the relaxation time. However, this 


modification requires evaluation of the post-shock levels of T, and T.. Interpretation of 
these quantities in a shock capturing solution is somewhat ambiguous [92). Improper 
definition can lead to instabilities in some circumstances, consequently no attempt 


was 


made to include this modification. 


Figures 7.4 and 7.5 present the effect of the different vibrational-translational en- 
ergy exchange cross sections in Eqs. 7.7-7.9 on the temperature and radiative emission 
profiles predicted for Case A. For a smaller tr„ the relaxation time r. increases. In- 
creasing r. delays the relaxation of translational to vibrational energy. This increases 
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Figure 7.4. Effect of <x v Models on Temperature Profiles - 1631 sec 

Tt while reducing T v , as illustrated in Fig. 7.4. Although the standoff distance in- 
creases with increasing T u the reduction in radiation due to the lower T v in this 
model more than compensates. The lower T v reduces the radiation from all three 
mechanisms, as shown in Figure 7.5. The lowest radiation for all three FIRE cases is 
predicted with Eq. 7.9, and the highest with Eq. 7.8. Comparing the predicted radi- 
ation for the vl 73, v273 , and v3T3 models in Tables 7.1-7.3 shows that the influence 
of o v is largest in Case A, the most nonequilibrium case. 

Again all three cases are closest to FIRE II when the radiation prediction is lowest, 
with Eq. 7.9, though only Case A exhibits much sensitivity or discrimination among 
the models. The predictions in Case B remain high, while the effect on Case C. 
the near-equilibrium case, is within the data uncertainty as was observed for the 
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Figure 7.5. Effect of a v Models on Radiative Emission Profiles - 1631 sec 
dissociation temperature models. 

7.2.3 Energy Exchange in Dissociation 

When a diatomic molecule dissociates, the vibrational energy it contained is con- 
sumed by the higher ground states (heats of formation) of the constitutive atoms. 
The energy thus removed must be accounted for in the vibrational energy equation. 
The amount of energy lost is commonly assumed to be the average vibrational energy 
at the local conditions, TT V . Because dissociation from a higher vibrational state may 
be more probable (the concept of preferential dissociation) several other formulations 
have been proposed [69]. 

One model assumes that the vibrational energy loss in dissociation, AE V , is some 
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fraction of the dissociation energy of a molecule measured from its ground state, D: 

AE v = c,D (0<c, < 1) (7.10) 

Sharma et al. [88] have proposed using ci = 0.3. This model was considered, but 
since it predicts the non-physical result of a negative T v in the shock region it was 

not pursued [92]. 

Park [93] has proposed that 

AE v = D-kT t (7.11) 

so that the energy lost is the dissociation energy minus the average translational 
energy. This model arises from the assumption that dissociation only occurs from 
levels that are within the average collisional energy of the dissociation threshold 
(similar to the mechanism of bound-free radiation). Again, this A E v is too large 
(particularly when T t is low) and results in an unrealistically small T v behind the 

shock [92]. 

A third proposed model assumes that dissociation occurs from some vibrational 
state(s) above the average. This is expressed as: 

A E v = c 2 E~ v {c 2 > 1.) (7-12) 

The best value of c 2 to use in this empirical model is unknown, and may depend on 
]? v ( 0 r T v ). In this work c 2 = 2 has been selected to provide an initial assessment of 
the model. The suffix 2tvs is used in the run identifier for these cases. 
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Figure 7.6. Effect of A E v Models on Temperature Profiles - 1634 sec 

Figure 7.6 compares the temperature profiles for Eq. 7.12 with c 2 = 1 (the baseline 
model) and c 2 = 2 for Case B. The increased A E v in the 2evs model (c 2 = 2) reduces 
T v in the region behind the shock where dissociation occurs, and increases it deeper in 
the shock layer where recombination begins. It also results in a smaller shock standoff 
distance. Figure 7.7 presents the radiative emission for Case B and confirms that the 
lower T v near the shock results in reduced radiation, while the higher T v in the shock 
layer increases it. This is observed in the tables as a decrease in the molecular band 
radiation and a slight increase in the atomic contributions, resulting in a net decrease 
in the total radiation. (The lower peak in the dashed curve of Fig. 7.7 is the atomic 

radiation peak.) 

For Cases A and B, the decreased radiation resulting from the 2evs model bring* 
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Figure 7.7. Effect of A E v Models on Radiative Emission Profiles - 1634 sec 

the prediction closer to the flight measurement, but it remains high (Tables 7. 1-7.3). 
For Case C the effect is within the data uncertainty, but this is the only model where 
the prediction is lower than the flight measurement. Table 7.3 also includes a solution 
run with the model of Eq. 7.10 using c x = 0.3 (denoted 3dis) for comparison. A flag 
was introduced in LAURA to suppress negative temperatures. The result is not 
significantly different so it provides no incentive to pursue this model. 

As with the two previous energy exchange models, the sensitivity to this model 
decreases as the flowfield tends toward equilibrium. The decreasing sensitivity to all 
these models near equilibrium is expected, since they govern the exchange between 
energy modes in the two-temperature thermal nonequilibrium model. When dealing 
with flowfields near equilibrium, the choice of formulation for these models is not as 
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important. The remaining models shown in Tables 7. 1-7.3 will be discussed later. 

7.3 Spectral Radiation Comparison 

To obtain further information on the validity of the energy exchange models, the 
FIRE spectral radiometer measurements (vs. wavelength, A, in fim) are compared 
with those exchange models whose predictions are closest to the flight measurements. 
Only continuum radiation is shown for the predicted spectra, since the resolution 
of the flight spectral radiometer is not sufficient to record any line radiation. In 
evaluating the flight spectra reproduced here, it should be noted that only a few are 
available in the literature. Those that have been published are the ones that were 
judged by researchers at the time to contain the least amount of noise. This criterion 
had the effect of selecting spectra that correspond to the lower range of flight radiation 
measurements. It is not known whether the unpublished data which correspond to 
the higher range of measurements had increased intensity levels or broader spectra or 

both. 

Case A (t=1631 sec) 

The spectral distribution of radiation has been examined for Case A for the v373 
solution, which appears in Table 7.1 to best match the flight data. Figures 7.8 and 7.9 
compare the prediction at 1631 sec with a flight spectral radiometer scan at 1631.3 sec. 
the closest available, reproduced from Cauchon [77]. Though the prediction is much 
higher than the flight scan they are qualitatively similar. Both show radiation from 
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the N * first negative band, whose two major band heads can be distinguished, (The 
spike in the flight spectrum between 0.5 and 0.6 pm has been attributed to a spurious 
signal [77].) Integrating the LORAN predictions over this limited spectral range 
results in a value that is high but within the data scatter at this flight condition [16, 

Fig. 13]. 

Case B (t=1634 sec) 

Figures 7.10 and 7.11 compare the Sees result tor Case B, the closest to the flight 
data in Table 7.2, to a flight spectrum at 1634.43 sec [77], the closes, time available. 
The qualitative agreement is comparable to that shown for Case A above. Integrating 
the predicted radiation over the range of the spectral radiometer gives a result that 
again, while high, is within the scatter of the flight data given in Fig. 13 of Cauchon 
[76]. The curve labeled dk73 in this figure will be discussed in Sec. 7.4. 

Case C (t=1637.5 sec) 

Figures 7.12 and 7.13 compare the predicted v37S spectrum for Case C to the 
nearest available flight spectrum, at 1636.43 sec. As in Case A, the v373 shows 
the closest match to the High, data in Table 7.3. By this point in the trajectory, 
the radiometer windows are approaching their melting point. A detailed post-flight 
analysis determined, however, that the windows cause less than a 10 percent change 
in the measured radiation [76], The intensity spectrum predicted by the v373 is 
excessive, having far too much radiative energy in the Jiff first negative band. The 
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Figure 7.11. Measured Radiation Spectrum - 1634.43 sec 
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curve labeled vleq will be discussed below. 

7.4 Nonequilibrium Chemistry 

The set of chemical reaction rates used in a flowfield model is somewhat bet- 
ter defined than are the rates of energy relaxation. Uncertainties remain, however, 
particularly in the extrapolation of the rates to temperatures at which little or no 
experimental data is available. Two rate sets are in widespread use: that of Dunn k 
Kang [94], and that of Park [87]. The latter rate set has been updated several times 
[67, 95]. In addition, there are several options for obtaining rates for the backward 
reactions. Dunn k Kang originally included curve fits for the backward rates. It is 
generally agreed, however, that more accurate results for backward rates are obtained 
by computing them from the forward rates and associated equilibrium constants (at 
appropriate temperatures). These equilibrium constants are commonly obtained from 
curve fits. Recently Gupta [96] generated a new set of curve fits for these constants 
which differ noticeably from those proposed by Park. Mitcheltree [97] has used the 
LORAN radiation model in the LAURA code to study the influence of these various 
sets of chemical rate data for a 12 km/sec entry. A few results for the FIRE II cases 

are presented here. 

The Dunn k Kang chemical kinetics model, denoted ik, has been applied to each 
of the FIRE II cases. As shown in Table 7.2, it results in a closer match of the flight 
heating rate for Case B. The wall spectral intensity predicted for this model was 
therefore also included in Fig. 7.10. This spectrum is much closer to the FIRE II 
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Figure 7.12. Predicted Radiation Spectra - 1637.5 sec 



Figure 7.13. Measured Radiation Spectrum - 1636.43 sec 
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Figure 7.14. Effect of Chemical Kinetics Model on Temperature Profiles - 1634 sec 

measurements than any of those obtained for Case B with the Park kinetics [87]. The 
distribution and relative magnitudes of the several band heads in this figure agree 
quite well with the flight data, and the Dunn k Kang prediction exhibits quantitative 
as well as qualitative agreement with the measured spectrum. For Case A, Table 7.1 
shows that the Dunn & Kang model is not an improvement. For the near-equilibrium 
Case C, the Dunn k Kang result was obtained for a thermal equilibrium solution. 
This resulted in a significantly higher radiation prediction (Table 7.3). 

Referring again to the temperature profiles in Figs. 6. 3-6. 5, it is clear that the 
thermal relaxation is in two stages with different time constants. Right behind the 
shock Tt decreases rapidly as a result of dissociation, then starts to level off. This sug- 
gests the involvement of slow reactions such as ionization, which become increasing!} 
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important as the level of ionization rises. This twostage equilibration is much less 
apparent for the Dunn St Kang chemical kinetics model profiles shown in Fig. i.H 
for Case B. The improved agreement noted (or Case B only with the Dunn k ha g 
solution suggests, therefore, that this second equilibration process deserves further 

study. 

7.5 Nonequilibrium Temperature 


The most obvious model of importance to the radiation calculation is the noneqm- 
librium temperature model. The choice of the number of temperatures modeled and 
the equations used is a crucial one for radiation. In a completely nonequilibrium 
situation, separate temperatures would be required for each energy mode of each 
species. This is impractical, except perhaps in a Direct Simulation Monte Carlo 
(DSMC) solution. Park's two-temperature model, discussed above, seems to provide 
reasonable accuracy with the minimum of complexity. If three temperatures are to 
be used, there is as yet no consensus about the electron energy equation which should 
be employed. Carlson and Gaily [981 have studied this question for Martian return 
conditions. Candler [4] has studied the problem for general conditions. 

For Case C, which is close to thermal equilibrium, solutions were generated for 
both the Park and the Dunn k Kang chemical kinetics, using Park’s two temperature 
model or assuming thermal equilibrium (a single temperature T). This allows an 
assessment of the extent and impact of thermal nonequilibrium. The results are 
included in Table 7.3 (vie, and dke,, respectively). The vie, result is closer to the 
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Table 7.1. Effect of Energy Exchange Models for FIRE II at 1631 sec 


Run 

Flux. W/cmMO.31-5.4 eV) 

Identifier 

Atomic 

Cont. 

Mol. 

Band 

Atomic 

Line 

Total 

flight 

vllO 

.IE-4 

.54 

.4E-3 

.16 ± 20% 
.54 

vl55 

.7E-5 

2.0 

.5E-2 

2.0 

vl73 

.IE-4 

1.1 

.9E-3 

1.1 

v273 

.IE-4 

3.2 

.2E-1 

3.2 

v373 

.IE-4 

.41 

.5E-3 

.41 

dk73 

.2E-3 

1.7 

.5E-1 

1.7 

vl732evs 

.IE-4 

.72 

.5E-3 

.73 


Table 7.2. Effect of Energy Exchange Models for FIRE II at 1634 sec 


Run 

Flux, W/cm 2 (0.31-5.4 eV) 

Identifier 

Atomic 

Cont. 

Mol. 

Band 

Atomic 

Line 

Total 

flight 

vllO 

.2E-1 

17.7 

2.7 

8.2 ± 20% 
20.3 

vl55 

.4E-1 

19.1 

3.8 

22.9 

vl73 1 

.3E-1 

17.7 

3.3 

21.1 

v273 

.7E-1 

22.5 

5.2 

27.8 

v373 

.IE-1 

17.6 

2.2 

19.9 

dk73 

.3 

3.1 

7.9 

11.4 

vl732evs 

.1 

12.1 

6.0 

18.3 


flight result and is therefore shown on Fig. 7.12. The agreement between this predicted 
spectrum and the flight spectrum is remarkable. (Recall that the spectral radiometer 
did not measure below 0.3 pm, so there is nothing to compare with in this range.) 
Thus, examining the radiation spectra at this condition suggests that Case C is in 
fact in thermal equilibrium, or at least that the distribution of energy between modes 
is more nearly in equilibrium than was suggested by the other nonequilibrium CFD 


models. 
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Table 7.3. Effect of Energy Exchange Models for FIRE II at 1637.5 sec 
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8 Conclusions 

8.1 Accomplishments 


A new method, LORAN, was developed to predict gas radiation in conditions 
of thermochemical nonequilibrium. This method includes a moderately de.atled ra- 
diative spectrum employing the smeared band mode, for molecular radiation. It is 
intended to provide radiation predictions with an accuracy between that of detailed 
line by line models and what is obtainable from highly approximate but very fast 
step models. The optimization of the method to provide the maximum detail for the 
minimum of computational effort is described. Representative results comparing the 
method to available ground and flight data indicate that LORAN predicts radiation 
with accuracy similar to that of Park’s NEQAIR code [20] and the nonequilibrium 

method of Carlson et al. [35]. 

Two options for radiative transport were incorporated in LORAN. The first is the 
traditional tangent slab method, in which the radiating medium is assumed to be 
one-dimensional. A second transport method was developed without invoking this 
approximation by applying a modified differential approximation (MDA). Predicttons 
of the radiative heat flux to the wall from the two transport methods were compared 
and show qualitative agreement with fluxes from the MDA method about 20 to 25 
percent lower than from the tangent slab. The profiles suggest that the boundary 
conditions of the MDA method might be improved. The variation of the divergence 
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of the radiative heat flux, appearing in the energy equation to couple the radiation 
and flowfield properties, was also examined and found to be much smoother in the 
MDA method. In addition to providing more accurate coupled results, this feature 
of the MDA method has the potential to enhance the stability of coupled solutions. 

A sensitivity study was performed using LORAN and flow-held solutions from 
LAURA to assess the effect that various models used in nonequilibrium CFD codes 
have on radiation predictions. The results indicate that radiation is a very sensi- 
tive indicator of phenomena occurring in nonequilibrium flowfields. This fact may 
be used in conjunction with currently available and future flight and ground-based 
data to improve the modeling of such things as exchange between energy modes in 
thermochemical nonequilibrium, and chemical reaction rates. 

8.2 Future Work 

Opportunities and needs for additional work have been identified in several areas. 

LORAN makes use of the NEQAIR method and data set for the computation 
of electronic energy level populations in nonequilibrium. This QSS algorithm incor- 
porates assumptions about the rate of change of energy level populations which are 
known to be violated in high gradient flow regions. The rate data which it uses are 
also known to contain uncertainties. Both of these areas of potential error should be 

investigated in more detail. 

Many of the radiative properties used in LORAN are uncertain, so the radiative 
heating values it predicts can only be regarded as approximate. An effort needs to be 
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undertaken to reduce these uncertainties. To do so may require new measurements to 
be made beyond those already available in the literature. In addition, corrections to 
the bound-free and free-free atomic continuum absorption coefficients accounting for 
deviations from a hydrogenic atomic structure could be included, though an initial 
inquiry suggests that these corrections have little effect. 

The multi-dimensional transport algorithm presented here is only a first iteration. 
Several aspects of this algorithm will be developed further in the future. These include 
the application of the method to cases when the absorption coefficient corrected for 
induced emission is nonpositive, the optimum selection of the underrelaxation param- 
eter for each spectral point, and the identification of improved boundary conditions. 

Additional work is needed to make LORAN a practical tool for computing nonequi- 
librium flowfields with coupled radiation, especially when coupling effects are signif- 
icant and many iterations are required to converge the coupled solution. Work is 
continuing to further reduce the time required to compute radiation by applying 
vectorization and other efficient strategies. 
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A Gaussian Units 

Much of the work in radiation has historically been done in the Gaussian system of 
units. This system introduces the particular oddity of measuring the electron charge 
in statcoulombs. As this particular unit of measure is uncommon, the following 
information is provided for the reader. 

In the Gaussian system of units: 


c=2.9979xlO l ° 

cm/sec 

speed of light in vacuum 

e=4.80286xl0 -1 ° 

statcoul 

electron charge 

h=6.6262xl0 -27 

erg-sec 

Planck’s constant 

k=1.3807xlCT 16 

erg/K 

Boltzmann’s constant 

m e =9.1095xl0 -28 

g 

electron mass. 


The following relations hold among several of these quantities: 
The Bohr radius, do, is given by 

an = — — = 0.52918xl0 -8 cm 

0 47r 2 m e e 2 

The fine structure constant, a, a dimensionless quantity, is 

9 tt6^ o 

a = — — = 7.292xl0 -3 


(A.l) 


(A. 2) 
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B Finite Volume Formulation of 
Radiative Transport 

B.l Finite Volume Development 

The following development is for the modified differential approximation for ra- 
diative transport as reported by Modest [271- The method has been reexpressed for 
a non-gray, non-scattering and nonequilibrium medium. 

The governing equations for the medium flux and incident intensity are: 

= 4 ^;' - ^ B-1 ^ 

SjG v = — ^ 

In what follows the R subscript on the radiative heat flux will be dropped for clarity. 
The first of these equations is a scalar equation, while the second is a vector equation. 
There are accordingly four unknowns for each spectral frequency chosen: G„, 9... 9... 
and . The only quantity required for coupling to a Computational Fluid Dynamics 
(CFD) solution, however, is the divergence of the radiative flux, V' 9R- Thls 
is given as a function of G by integrating Eq. B.l over frequency, so that a solution 

for the single unknown G„ at each spectral frequency is sufficient. 

To obtain a form of these equations suitable for a flnite-volume solution, first 


-i-VG. = -39. 


rearrange Eq. B.‘2: 


(B.3) 
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Now take the dot product of this equation to allow the substitution of Eq. B.l. 



= -3y • & = -12 jt jI + 3 k' u G v 


(B.4) 


So now rearranging, 


v • ^vG^ = 3 ' C ‘ /G " ~ 127r - ? ' 


(B.5) 


In the finite volume approach, this is now to be integrated over the volume of a single 


grid cell, V. 

j I j dV = / / j v (3<G^ - 12tt jl)dV « V(3k ' v G v - I2ir ;*) 

(B.6) 

Then by application of Gauss’ theorem (the divergence theorem), the volume integral 
on the left-hand-side can be transformed to a surface integral. 


J J ^-^yG* • nj da = V {Zk' v G u — 127r;*) (B.7) 

The integral is performed by assuming that the absorption coefficient < and the 
gradient of the incident intensity V G * are constant on each cell face. The integral 
then becomes a summation over the six faces of the cell. The notation used for the 
cell geometry is that of Gnoffo [99], where /, J , and K denote cell centers, and i, j, 
and k denote cell faces. Then 


(-tVGv) -h.+ia.+i 
[\< A+i 

• ^j+i a j+i 

\< ) J+l 



n,a. 


J J,K 



Tljdj 


I,K 
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+ 



hk+l a k+l 



tlkO-k 


J I.J 


= V [Zk'^Gu — l'-Trj')/ j k 


(B.8) 


For the nongray gas found in a shock layer, the radiative properties and j u 
vary over orders of magnitude at various locations and frequencies. To minimize the 
numerical difficulties that this variation can introduce, some kind of normalization is 


desired. Define 

•? _ Ecells jt (B.9) 

^cells 

and 

_ _ EceUs K 'v (B.10) 

^ ^cells 

Both these quantities are constant for each frequency considered, and provide a mea- 
sure of the approximate magnitude of the emission and absorption at that frequency. 
Figure B.l is an example of the variation of these average coefficients for the Mars 
return test case presented in Ch. 6 and demonstrates the wide variation in radiation 
properties for various frequencies in a single flowfield. 

Proceeding with the finite volume development, now divide Eq. B.8 by J„ on both 
sides, and also multiply terms containing G u by 1 in the form of TS w /k„. Define a new 
variable r„ = G v lQ y K u ). Then Eq. B.8 becomes: 


— vr* 

k k' 


hi+xdi+i — , vr, ) * rc i a i 


*+i 


J,K 


+ 


— vr. 




Aj +l aj+i - ( vr A ■ rijdj 

Ilk 
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Figure B 


. 1 . Variation of 7c„ and j u for a Mars Return Flowfield 


+ 


^fvr. 


nk+\ a k+\ 


k + 1 


-(!*•). 




i' 




= V ( 37c, - 12tt^ 


[B.n; 


J v / I,J,K 


This formulation reduces the variation of the unknown, T„ and of the coefficients of 

the equation, making the numerical solution easier. 

Now consider the j and j + 1 cell faces in an axisymmetric flow. The gradient 
of r, has no component in the circumferential direction, while the surface normal n 
on these faces is entirely in the circumferential direction. Therefore the dot product 
of these two quantities is zero and the terms on these faces drop out. This is a 
mathematical expression of the fact that there is no flux through these cell faces in 
an axisymmetric flow. The j and J subscripts then become superfluous, and will be 


omitted in the remainder of this development. 
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Since the flowfield grid is not orthogonal, it is necessary to use a generalized trans- 
formation to obtain finite difference expressions for the gradients in Eq. B.ll. The 
elements of such a transformation are repeated below from Anderson et al. [100, Sec. 
5-6.2], specialized to the case of an axisymmetric flowfield. Second order accurate cen- 
tral difference expressions are used to obtain the necessary partial derivatives, where 
£ r,, and C are the coordinates of the orthogonal and uniformly spaced computational 
grid corresponding to the i, j, and k grid directions, respectively. 


where 


x;+i ,jt - a?i-i,fc 
“ 2 

(B.1'2) 

X-q 0 

(B.13) 

Xi,k+ 1 “ X i.k- 1 

x C — 2 

(B.14) 

o 

II 

uy 

5* 

(B.15) 

y.,i+i,fc - — x, k sin 0 

(B.16) 

9 = 2. * arctan -| ce n CO mer 

X 

(B.17) 

O 

II 

(B. IS) 

Zi+l.k — -1—1,* 

2 

(B.19) 

^ = o 

(B.20) 


(B.21) 


= 2 
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At the boundaries firs, order backward differences are applied, excep. for the axis 
Of symmetry boundary. At this boundary the coordinates of a “ghost” cell can be 
obtained through reflection. Second order accurate central differencing can therefore 
be applied at this boundary. The Jacobian of the transformation can then be obtained 


from: 


1 

(B.22) 

J 

— x (,Vn z i 


while the metrics of the transformation are given by: 


6 = hvH 

(B.23) 

o 

II 

w 

(B.24) 

i z = -JxcVr, 

(B.25) 

H 

II 

O 

(B.26) 

T) y = -X C Z € ) 

(B.27) 

N 

II 

o 

(B.2S) 

N 

cr 

1 

II 

O 

(B.29) 

Cv = 0 

(B.30) 

C* = Jx&rt 

(B.31) 

The first derivatives which appear in the gradient terms are 

then obtained from 

II 

+ 

J"', 

(B.32 
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d 


d 


dy ' h dy 
d . d . d 
d~z = ^d^ + ^d c 


(B.33) 

(B.34) 


so the gradient can be expanded as: 


- / dT u . <3IV\ * dT u - ( dT y . dTA ? (B 35’ 

vr* = ( £x^- + ) * + ny—J + + <>•— 1 k 


dr, J 1 V sz K d < 

The y-derivative can be omitted since there is no circumferential variation of IV in 
the case of axisymmetric flow. Collecting terms then yields 


vr. = (t! + a) f + (oi + c-t) f + vcf (b. 36, 


Now note that the combination ha is just the directed area of the cell face, or a. 
The finite volume expression, Eq. B.ll, can now be reformulated as. 


(vf. ■+>.*• +V0+1.A— ) 

«!. . . V a K i +\ ,k 


>i+l, A' 


Kv (vtijc^r + VC..A' 


,+1,/v , 


A 


^/.fc+l 


1C!.... V at,* ' — dc iiK 

dT u _ . 

,*+i 


• a., K 


Vi.K 

f - e dl„ 
V?/,H 1 


+ vC/.h-i -ttt )‘ a ^+i 

^ /.A 


,fc+l> 


dC /, 

~k u ( - e dT „ - 5T 

-7- + VC/.* 

v 




9IV \ 

<9( /.J 


• a/,fc 


= V/,a- 3k„icX - 12jt^ 

V A/ i,K 


( B . 3 7 ) 


The radiative properties < and j* are known at cell centers The required 

values of k[ at the cell faces can be obtained by averaging the values at the two 
nearest cell centers. Then for instance 


1 >/ * - + V 


K U+l.K “ V^Z+l.K ^I.K, 


(B.38) 
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This average has been found to produce better stability than the alternate form with 
+K > ) because of the extreme variation of the absorption coefficient near 

the freestream boundary in particular. The first derivatives of T. in computational 
space at the cell faces are obtained by second-order accurate central differences, as 

follows: 

dY u 


gr, 

d( t+i, a' ^ 


ai u _ /r - r ) 

-TT- — U^+l.A- Vv 1,k) 

di i+i,K 

1 K±l ~ EfZ£±l _ ^ l//,A '- 1 


(B.39) 

(B.40) 


The next tw 


m differences are the same but with a shift in the i and I indices: 


d 4r =(r, f , K --i W) 

di i,K 

dT „ 1 ( Tuui+iJzSjZJSzl. | r ^-' JLtl ~ r ‘ /f ~ 1 

5Cijir“ 2 \ 2 2 ' 


(B.411 


(B.42) 


Along the normal cell faces: 

d T u 


1 (Tvj+i k + i ~ | ^ni.K (B.43) 


i,k+i ^ 


— tr — r ) 

d( I Ml 


(B.44) 


And again the next two 
indices: 


differences are the same except for a shift in the k and A 


dr 

dt 


, _ l / r, f+1K t ^±le=i _ 

; i,k - V 


dr v 
dc, u 


— (Ti'/.K r't'/.K-l ) 


(B.45) 
( B .46) 
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All these expressions can now 


i be substituted back into Eq. B.37, to obtain: 


i'l + l ,K 


vij< * 




1 -r, f+1 ,K-, + ^k+i — tiLhsi 


y -~ — 

* CE t + 1 ,/v 


o 


p- [v6.A'(r^;.K _ 


v t,I< 


1 (T 


E2=L £=1 




+V4i,K'9 i 2 

i / r„ +1 K + l 


+ 


2 




t'J.fcH- 1 


VCu+lTj 


- r^.j.A+1 , jVtiw 

O 


+ VC/4+l( r ‘'I,A'+1 - r ^/.A-)] ■ a/ ' fe+l 


i -L Ei2±i^=i — ^ l ' J - 1 — 


t'/.k 


VO,fe 9 

— T )] • «/ fc = V/,k f3«u K i/Fi/ ~ 127Ty- j 
,k- 1 *'/>•-» ; J J ’ \ Jv/I,K 


Carrying out the dot products and collecting terms results in: 


hr 


t'J + l.fC+l 


+ + -MV,.,*-. + 


+*r^„ + - - 12 ’ rV «X (B ' 48) 


(B.47 


where 


VCi+i.K’ ' a.+i.ft' + V6,fc+i • a ^+i 

4 K ' . . „ "M+i 


1 


J't+l.K 


(B.49) 


— 1 K u - £ 1 H. jj. t—y£ r r ■ fl J fe (B.50) 

d 2 = -A-V«. + r.a- • S. + .,K + J— Vfnw •«'*« -K,, , V{ “ 

^ ^ y/.k+l 


I'l + l.K 


d3 = -i^VC +1 .K-^.-^V6,- 3 u 


I'l+hK 


(B.51) 
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1 

d * = 1^- 


'‘'i+l./C 


V0+1.A- ■ feur - J^-VOA • ** + ' toM ‘ B i2 » 


d 5 = — -^— v£;+i,a- • «<+i.a' - zr~^^' K ' a ‘ K 

< + l.K 


— VCu+i • «M+1 - -T^VCu • - 3VU-K.K 

<x,* +1 K *"* 


"/.A' 


* = ■ -- + ■ *•* + t;* <M • Sm (B,54) 


1 hy 


i'l+l.A' 


d 7 = -7-t^VC.,k- • a., A' - ’ VO,fc+i ' “/.*+! 

4 V/c 


I'/.k + l 


d» = 4^-V{i,A • Si.K - J^VO.W • + i77 Vf, '‘ ' Jw (B ' o6) 

/C.. „ ** * */,*+! * 7 ’* 




(B.ol 


d 9 = j-T^VCkA • ai, A + ‘ Sl ' k 

4 <,k- 4#£ ^ * 

The dot products of the metrics with the cell face directed areas which appear 

in these coefficients are determined by the grid geometry. These can be written as 


follows: [99] 


V&',A • a., A = 


• u;,A = 


+ 


V, K^i-hK + flj.fr) + Vf-uc( g <.A + dj+ij<l 
4V7 ,a V/-i,a 

l^Vf, K + + fl ^..K a x,K- + ^fj£> 

4 V/,a'V/--i,a 

V7-l,A(Qr,,K Q ^. K + a --..A a -. K + a «.-H -A ag CA- + 

4V7 i r-V;_i,a 


(B.58) 


[ V / /.A(a..A-l + Q.-H.A-l) + V7.A-l( a »,A + filial | , q / fc - 

Vs/,/: • O/.fc - 4V7.A-V7.A-1 J 

V7.A(flx i .K--i a r/,* + + a * .-n.K-i a *;.> ± 

4V7,a'V/,a-i 

V7 K-lK t + a *,.K a H.k t fl ^».A- fl */ ■* + 

4V7,a17,a-i 


(B.59) 
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• ai.K 


V / , K (g/- 1 ,fc + ai-i.k+ 1) + Vi-iM* + £[±tll 
VC./v ’ a i,A — 4 V[ t K V 

V, h -{ a X ' .Qx, A~ + !i + ^Ilz. L*±lf^ 

4V/,a ■ V/-i,a 

V/-1 ,K( a Jj ), q r, A- + a zi.k a -,.K + a rj,M-l Q j - ./v + 

4V/,k^/-i.k 


(B.60) 


• a/.jfc = 


V[,l\( a I,k-l ~6 Qf.fc) d* W.K — 

VC/.fc ■ — 4 V/,kV;,K-1 

4V / /,k V/.K-l 

V7.y-i(qg/.* a g^ + q «/.* a »/-* + + Q «/.*+» a */-*I 

4V/,/ v V/ i a'-i 


+ 


( B .61 ) 


Equation B.48 is a matrix equation for the unknowns r„ of the general form 
Ax = E. The matrix A and the vector 6 are functions only of the geometry and the 
absorption and emission coefficients. The numerical solution of this type of equation 
can be obtained by a number of methods, including direct inversion of the matrix, and 
various relaxation methods such as point Jacobi, point Gauss-Seidel, line Jacobi or 
line Gauss-Seidel algorithms. The line Gauss-Seidel algorithm has been selected and 
applied along each normal grid line to capture the dominant gradients in the radiative 
properties. In this approach, all the T, terms along the I grid line are at the n + 1 
time level, while all / + 1 terms are at the n time level. Terms with I - 1 are known at 
time level n+1 from the solution of the previous line. The d coefficients, as mentioned 
above, are functions only of the geometry and the absorption coefficient. Since both 
these quantities are constant for a given transport solution, the d coefficients need 
only be calculated once. 

The line Gauss-Seidel algorithm results in a tridiagonal matrix equation to be 
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solved for each normal grid line. The solution of such an equation can easily be 
obtained using the Thomas algorithm [100] with specification of appropriate boundary 
conditions as discussed below. To ensure convergence for this problem, which is a 
nonlinear, elliptic equation, underrelaxation is recommended. In fact it is found to 
be necessary m this application. It is introduced by defining a corrected value for the 

update of r„ & s 


r~J' = r" ,. + r(r;«. - r;, K ) = d - dr;,,,. + -t 


>n+l 

t'J.K 


(B.62) 


where r is the relaxation parameter and is less than one for underrelaxation. The 
expression for r~J is obtained by solving Eq. B.48. Incorporating this expression in 
Eq. B.62 and rearranging to the line Gauss-Seidel form results in 

j 1 lp+1 + P +1 ' + = (1 - r ) r “, v _ X + 1. 

+d 3 r: /+1 K _, + dr r£\ iK+1 + d 8 r£\' K + *•-» + 127rV/ ' A '7T ) B-63) 

The selection of r is discussed in Sec. B.3 below. 

B.2 Boundary Conditions 

Figure B.2 shows the essential features of the radiation grid. Properties at any 
interior point (denoted by filled circles) can be obtained by the solution of Eq. B.63 
above. The unfilled circles at the left represent -ghost” points for the symmetry 
boundary condition at the axis of symmetry. Symmetry boundary conditions also 
exist in and out of the page on the -pie slice” sides (see Fig. B.3). bu, these are 



14S 


implicit in the axisymmetric formulation. The triangles along the upper edge of the 
grid represent the outermost grid cell which forms the freestream boundai\. The 
downward-pointing triangles along the lower edge at right are "ghost points for the 
outflow boundary. Finally, the barely visible unfilled squares along the lower edge 
of the grid represent the center of the surface elements. Each of these boundaries is 
discussed separately below. 



Figure B.2. Essential Features of Axisymmetric Radiation Grid 

B.2.1 Axis of Symmetry Boundary 

This is a symmetry boundary, at which reflection boundary conditions should 
apply. Therefore for the unfilled circles in Fig. B.2 along this boundary: 


r^ghost, K) = r„(i, a ) 


(B.64) 
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Figure B.3. Top View of Axisymmetric Radiation Grid 
This boundary is a singularity in the computational grid, however, where one 


surface of a grid cell has zero area 


on 


the second line disappears. The two remaining 


. Referring back to Eq. B.37, this means the term 
^-derivative terms at the I grid 


line are replaced by first order backward differences to avoid differencing across 
singularity[92]. Then 

5r i/ ^ ( v r 

u? /.fc+i 


the 


dt 


+ r„ /+1 . K e„ / iK ) 

(B.65) 

'*//+!. K-l _ E^,K-l) 

(B.66) 


/ , Ar 


Substituting, Eq. B.47 becomes 


■ — [y$i+i,A'(E/ + i,K Eu/j,' ) 


I'. + l.K 


- i /r „ f+liJf+1 Et / J+1 ,/v-i 

+VU1A9 1 2 




ni+i,K 
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+ - 


"/,) t+i 


+vCa(r,, K - r„ fiJf _,)] V 3 '^ 1 '" " l H, 


va* + 4(^K> l -W>. + r ^*.« r ^' 

+ V(f,Hl( r ‘'/,K + i ' a ^fc+l 

Vs/A-^ (r^r+i.K _ + ^/-h.k-i ri/f K -‘) 


(B.6T 


l y K 


Carrying out the dot products and collecting terms 
be put in the form of Eq. B.48, with 


as 


before allows this equation to 


1 K “* -» 

1 Ki/ -/• Z. , . - A - Vfj k+l ' a /,fc+l 

4 = 7- VUl.K • a l+ i,K + 2 , w '* + 

4 K 


(B.68) 


^» + l,K 


I'J.fc + l 


do — 


1 7c„ - f - _ I JlZ-vilk ■ aik (B.69) 

^ - c „ . 5... *- + — V6,fc + 1 • a ^- fc + 1 9 *' w ’ 

Vsi+t.K a.+i,A ^ 9 / ^ 

“ I'J.k+l 


^t + 1 ,K 


,3 = -li-vfaA • w - ' * M 




1 ^ r r., 
.1 

l »I,k 


(B.70) 


d.A — , 


1 Ku 


1 ic„ 


4 ic 


■VCi+l.^ ■ a * +1 ’ K 2/c' 


VU,fc+i • a ^+! 


I't + l ,K 


‘'J.lc+l 
K 


+ / — VO.fc+1 ■ a/ ' fc + 1 




(B.71) 


K - _ I_ Kl/ 

4 = - 7 — Vfr+1.* ■ a '+ 1 ' A ‘ _ 2 ic' 

fc*' L 


-V^,fc + 1 ■ a/ ' fc + 1 


^t+l,K 


^/,/c+l 




1 «u 


/i p 


"Uc + l 


•VC/,fc+i ’ g/.m-i + 9 it* 1 


V6,fc • «/.* “ 77 ’ a ^. fc 

— ‘XVi,K^v^vi,k 


(B .72 ) 


ft p> ~* . 


4 = -I-^-vC+n* • + 2^" V °’ fe ' GU + a/ " 

4 K 


t'l + l.K 


*/,fc 


(B.73) 
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d- = 0 


(B.T4) 


d s = 0 

^9 = 0 


(B.75) 


(B.76) 


With the appropriate values of d, along this boundary, the 
the same method as the interior points. 


be solved with 


B.2.2 Freestream Boundary 

The freestream “boundary” is transparent to radiation, and only outgoing radia 
,io„ is considered. The freestream is assumed to be non-emitting. With these assump- 
tions, this boundary can be approximated as a cold wall with complete absorption 
( £ „ = 1.0). The modified differential approximation is essentially a P-1 method, for 
which cold wall boundary conditions can be obtained analytically. Modest suggests 
the use of Marshak’s boundary condition, which is 


_ 2 [A-l] g Vn + Gv = 0 (B ’ 77) 

where *, is the spectral surface emissivity, and A is the outward unit normal vector 
describing this boundary surface. Replacing q„ in this equation using Eq. B.2 leads 
to an equation containing only the one variable G„: 

2 [i_lli% .ft + Gc = 0 < BJS) 

ISu J 

Then dividing by Q,H.) to obtain an equation for T, finally gives 

’IHft 


■ h + r„ = o 


(B.79) 
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- . dr, - . . dT„ 

v?n aj" + v; '" ac 


+ r„ = o 


B.SO) 


Expanding the gradient of IT gives: 

3< Vi„ 7 

This can be differenced by expressing the derivatives of IT using a central diffeience 
for the ^-direction and a first-order backward difference for the C -direction. Denoting 
the freestream boundary grid cell by the subscript F, the result is 

|V£/.F • u/.F- r ^ l F o + VC/.f ‘ hl ' F ( r,// F " 


3k 


,F 


+r„, F =o (B.si) 


where the emissivity for the freestream boundary was set to 1.0 as mentioned above. 
Collecting terms and arranging them in the line Gauss-Seidel form, with the addition 
of underrelaxation, gives 

-2t-vG,f ' nr.rC+J., + ( 2 VG.r • W + K,,) r -5 = 

(2VO.F ■ n,. F + 3<„) (1 - r)C, - C-',., > < B ' S2) 

The do. product terms calculated in Eqs. B.59 and B.61 are a. cell faces. To obtain 
the values at the cell center required in the above expression, an average at the ttvo 
opposing faces of the cell is taken. Then 

V f/.z+i ' n/j+1 + Vfr./ ' ± Vfr+W +i ' TTtH±l + VO+W ' Hlj±L 

V^.F ‘ fif.F = 4 

( B . S3 ) 

V C/.f+i • n u+ 1 + VOj • »/./ + VC/+ 1./+1 • ni+u+i + VC/+U • ? DtiI 
VC/.F • fi/.F = ' 4 


(B.84) 



where the subscript / + 1 denotes the freestream cell tace. 

For the freestream boundary at the axis of symmetry, the reflection boundan 
condition Eq. B.64 for r„,_ liF is used. Using a first order backward difference at this 
point was found to lead to nonphysical results. 

B.2.3 Outflow Boundary 

Along the outflow “boundary”, r„ /+1 K in Eq. B.47 represents the "ghost point 
and must be approximated. This is an arbitrary computational boundary through 
which radiating gas and radiation (inwardly and outwardly directed) both pass. The 
flowfield boundary condition used in LAURA is a zeroth order extrapolation (i.e.. 
the derivatives of the flow variables are assumed to be zero across this boundary). 
This boundary condition is known to provide better stability than other possible 
formulations, and so can also be adopted for the radiation method. Therefore, 

r„(/ + i,/0 *r„(/,A') ( B - 85 ) 


One alternative possibility is to assume a constant slope at this boundary. To perfoim 
an extrapolation of the positive-definite variable r„, a logarithmic extrapolation is 

indicated 1 


r. ( / + i.A-)« r 


(B.86) 


This boundary condition was found to work well in some cases and to cause T^(/. A ) 
to blow up in the shoulder region in others. This undesirable behavior results from 


Suggested by Gnoffo 
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inaccuracies and lack of grid resolution near the shoulder of the body in the test 
cases. Logic was implemented to apply the logarithmic boundary condition except m 
cases where it would result in an increasing value of T „{I,K) around the shoulder. 
Those instances are treated with the zeroth order boundary condition to prevent the 
solution from blowing up. Even when the solution does blow up at the shoulder, much 
of the rest of the flowfield remains quite unaffected. This suggests that the solution 

method is robust. 


B.2.4 Wall Boundary 

For the medium only intensity, this boundary appears as a cold wall. The Marshak 
boundary condition (Eq. B.79) can be differenced at this boundary by recognizing that 
the grid lines are normal to the body surface (Fig. B.2 is not to scale). The term 
VT, • n is then simply the gradient of T„ along the grid lines. This difference is best 
expressed in physical coordinates with a first order forward scheme as 

IY; i - rVu, 3/c„ ri ( £ *i \ r — Q (B.87) 

Ar 2 \2 — £ U [) ^ 

where the subscript w denotes a “ghost” point on the surface of the body, and 

Ar = y/{x{L 1) - x{I, u)) 2 + u; )) 2 (B - SS) 


In the Gauss-Seidel format, with the incorporation of underrelaxation, this becomes 


3/c' 


1 

Ar 


-i// 


9 


9 - r 


l 'J . 


r n + 1 ' - r — r n+1 = (1 — r) 

L vi. w Ar* "M V ' 


Ar 


Ar 


3k 


u i,i 




9 - r . 


'VI 


r n 

V[, W 


(B.89) 
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B.2.5 Radiative Heat Flux to Wall 


The radiative flux directed to the wall due to emission in the medium can then 


be found from solving Eq. B.2. 


ok' 


(B.90) 


In this equation, is obtained from the solution of the finite volume problem, T^, 
multiplied by the normalizing factor (] y K v ). To obtain the flux directed at the wall, 
the dot product of q v with the wall-directed surface normal is required. Then 


Quiv — 9^ ' <j^/ ^jGi, ’ 


fB.91] 


Again recalling that the grid lines are normal to the body surface, the gradient of G„ 
can be expressed with a one-sided difference along the normal grid lines. As above, 
this is best done in physical coordinates. 

1 G uj 2 G UI1 




(B.92) 


3k', Ar 

Alternately, a second order forward difference could be used to represent the gradient. 
The total heat flux to the wall at each grid location is then obtained by a numerical 
integration over the spectral result above. This closes the numerical problem. 


B.3 Numerics and Convergence 

B.3.1 Convergence Criterion 

The numerical solution of this set of difference equations requires the selection of 
an appropriate criterion to test for convergence. Because the radiative properties vary 
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over many orders of magnitude, and because the grid introduces numerical errors, the 
usual definition of the I 2 norm for this purpose is not satisfactory. Instead, a local 
error function is defined: 


ir;,y - 

mean(r") 


(B.93) 


where the mean of r* is obtained as the average of the minimum and maximum values 
of r n for all / and K on the radiation grid. This definition reduces the value of 

U I , K 

l l K in regions where V n UIK is small, relative to the error that would be obtained from 
the I 2 norm. The result is that convergence is reached faster. Though convergence 
may not be entirely complete at those locations where K . is small, these locations 
contribute little to the coupling with the flowfield. Therefore, the lack of complete 
convergence at these locations is deemed acceptable to reduce the computation time. 

The actual convergence is determined by the average value of F/,a' over the radia- 
tion grid. When this average is less than some specified value, convergence is assumed 
to have been reached for that particular frequency i/. Each frequency is converged 
individually, since the rate of convergence depends on the magnitude of the optical 
depth and varies considerably with frequency. A typical convergence history for a 
single frequency is shown in Fig. B.4. Convergence may be faster or slower for the in- 
dividual frequencies, depending at least partly on the magnitude of the optical depth 


at that frequency. 
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Figure B.4. Typical Convergence History for MDA Solution 

B.3.2 Selection of Relaxation Parameter 

The relaxation parameter required to obtain a stable solution varies with the mag- 
nitude of the optical depth of the flowfield. Very fast convergence can be obtained for 
optically thick spectral regions using r=l, meaning no underrelaxation. For optically 
thin regions on the other hand, the solution with r=l is unstable. The value r=0.o 
has been selected as a good compromise. In future, an algorithm might be developed 
to vary the underrelaxation parameter r as a function of the radiation properties in 
order to further accelerate the convergence of the complete transport solution. 


